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ABSTRACT 

We develop an algorithm for setting up initial Gaussian random density and velocity fields 
containing one or more peaks or dips, in an arbitrary cosmological scenario. The intention is 
to generate appropriate initial conditions for cosmological N-body simulations that focus on the 
evolution of the progenitors of the present-day galaxies and clusters. The procedure is an application 
of the direct and accurate prescription of Hoffman & Ribak (1991) for generating constrained 
random fields. 

For each peak a total of 21 physical characteristics can be specified, including its scale, posi- 
tion, density Hessian, velocity, and velocity gradient. The velocity (or, equivalently, gravity) field 
constrants are based on a generalization of the formalism developed by Bardeen et al. (1986). The 
resulting density field is sculpted such that it induces the desired amount of net gravitational and 
tidal forces. 

We provide a detailed mathematical presentation of the formalism. Afterwards we provide 
analytical estimates of the likelihood of the imposed constraints. Amongst others, it is shown that 
the tidal field has a strong tendency to align itself along the principal axes of the mass tensor. 
The method is illustrated by means of some concrete examples. In addition to the illustration 
of constraint-field correlation functions and how they add up to the mean fields, followed by il- 
lustrations of the variance characteristics of field realizations, we concentrate in particular on the 
consequences of imposing gravitational field constraints (or, equivalent in the linear regime for 
growing mode fluctuations, peculiar velocity field constraints). 

Subject headings: Cosmology : theory - Galaxies: clustering - large-scale structure of the Universe 
- Methods: numerical 

1. Introduction 

In the standard scenario of structure formation galaxies and the large-scale structure form through 
the growth of primordial density perturbations. These perturbations take the form of a homoge- 
neous and isotropic random process. In most cases these cosmological density fields are assumed 
to be Gaussian random fields. 

In these density fields the regions around local maxima and minima are of particular interest 
during the evolution of the perturbation field. The first collapsed structures form generally near 
(but are not coincident with, Bertschinger & Jain 1994) density peaks, making density maxima 
the progenitors of objects like galaxies and clusters. On the other hand, the minima will be the 
centres of expanding voids. The properties of peaks in Gaussian random fields have been described 
extensively in the literature. In order to identify an object of a certain size and mass in a Gaussian 
random field one discards smaller scale objects from consideration. This is achieved by filtering the 
density field on an appropriate scale to reflect the linear evolution of the proto-objects. While in 
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some scenarios the filter function is a consequence of a simple phenomenon (e.g. free-streaming of 
neutrinos in a Hot Dark Matter scenario) in other cases one is forced to invoke an artificial filter 
to approximate the complicated processes of hierarchical merging (e.g. in the Cold Dark Matter 
scenario) . 

A description of the properties of these filtered fields was given by Doroshkevich (1970), Peacock 
and Heavens (1985), and Bardeen et al. (1986; hereafter BBKS). Beside global parameters such 
as the number density and spatial correlations of peaks found in these filtered fields they also 
derived the distribution of their height, shape and orientation. Furthermore, BBKS derived the 
mean and variance of the density profiles around peaks. As soon as these structures enter the 
nonlinear regime the coupling of modes breaks down the above approach of filtering. To investigate 
the further evolution one is therefore forced to resort to N-body simulations. However, in order 
to follow the evolution of a particular object one needs to be able to start off with a primordial 
density field containing such an object. 

Unfortunately, the methods of BBKS apply only to point processes and cannot be used to 
construct an actual sample of a density profile around a peak with predetermined parameters such 
as peak height, shape and orientation. The usual approach is therefore to generate an unconstrained 
realization of a Gaussian field and then to search for peaks or regions that satisfy the desired 
constraints. In many instances this is an inefficient approach. For example, giant clusters or voids 
will be so rare that either many samples have to be generated or that a large box needs to be used to 
ensure that the object is indeed present in the simulation volume. The latter will yield a severely 
degraded resolution which conflicts with the desire to describe these objects in as much detail 
as possible. Similar considerations apply when many properties need to be specified to obtain the 
desired object, even while the corresponding additional constraints do not represent unlikely values. 
By being able to specify beforehand some of the properties and to ensure the presence of such a 
peak or region in the simulation volume the required effort will be minimized. Simultaneously, the 
resolution will be maximized. Potentially the most important advantage of this approach is that the 
influence of several physical quantities on the evolution of structures can be studied systematically 
by generating realizations wherein one or more constraints have various values. 

The fundamental theory of these constrained random fields was set forth by Bertschinger (1987). 
He generalized the treatment used by BBKS to give a full statistical description of a Gaussian 
random field subjected to constraints. Based on these principles he presented a method to correctly 
sample the probability distribution of the density field subject to linear constraints. This method, 
however, is rather elaborate and inefficient in its implementation, involving a simulated annealing 
technique. Although it is useful for generating initial conditions subject to a few constraints (see 
e.g. Van de Weygaert & Van Kampen 1993), it quickly becomes prohibitively slow for more than 
two constraints. Looking for a more efficient procedure, Binney and Quinn (1991) showed that 
Bertschinger's problem simplifies considerably when the random field is expanded in spherical 
harmonics rather than in a plane wave basis. In the case of a localised set of constraints, such as 
the presence and shape of a peak at the centre of the box, the problem can then be solved exactly 
instead of iteratively. However, their algorithm is essentially restricted to the case of quite localised 
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constraints denned around an obvious centre of symmetry. 

The breakthrough in the construction of constrained random fields came with the publication 
by Hoffman & Ribak (1991, hereafter HR). They realised that for any constraint that is a linear 
functional of the field the problem can be solved exactly in an elegant and simple manner, without 
having to invoke complicated iterative techniques. Their method makes it possible to generate 
initial conditions for N-body simulations that obey a few hundred constraints, e.g. those imposed 
by the observable universe (see Ganon & Hoffman 1993). 

This paper contains a description of the fundamentals and implementation of a specific cosmo- 
logical application of the method proposed by Hoffman & Ribak (1991). This application consists 
of the generation of an initial density and velocity field containing one or more density peaks in 
a simulation box. Apart from being able to determine the location and the scale of the peak, we 
can specify the central density of the peak, as well as the compactness, shape and orientation of 
the density field in the immediate surroundings of the peak. In addition, the total matter distri- 
bution can be sculpted such that it subjects the peak to a desired amount of net gravitational and 
tidal forces. In practice, the computer algorithm generates samples of these constrained Gaussian 
random fields on a lattice, using Monte Carlo techniques. Nearly all relevant calculations are done 
in Fourier transform space. Some results of cosmological studies based on these constrained initial 
conditions are presented by Van Haarlem & Van de Weygaert (1993), Van de Weygaert & Babul 
(1994, 1995). 

In this paper, we start with some basic concepts of Gaussian random fields followed by a 
treatment of the fundamental theory of constrained Gaussian random fields in section 2. The 
Hoffman- Ribak method for the construction of constrained random fields is described in section 3, 
followed by a description of our Fourier space implementation. In section 4, we present our appli- 
cation of this method to the generation of peaks, deriving constraint kernels for the various peak 
quantities. In addition, we provide prescriptions for the probability of the imposed constraints. A 
realization of a random density field with a constrained peak is presented in section 5. Specifically, 
we will focus on the influence of imposing a peculiar acceleration and a tidal field. In section 6, we 
will conclude with a summary and short discussion. 

2. Fundamentals of constrained Gaussian random fields 

Although the paper by Hoffman and Ribak presents the essentials of the simple direct method to 
construct samples of constrained random fields, it does not provide its mathematical background. 
This can be obtained extending our earlier treatments (Bertschinger 1987, Van de Weygaert 1991). 
Therefore we will first summarize the necessary mathematical background in the notation employed 
by HR before we get to the presentation of their method. 

2.1 Gaussian random fields: basics 

Consider a homogeneous and isotropic random field /(x) with zero mean. The random field is 
defined by the set of Appoint joint probabilities, 

V N = P[/( Xl ), /(x 2 ), . . . , /( XJV )] d/(x 1 )d/(x 2 ) • • • d/( Xjv ), (1) 
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that the field / has values in the range f(x.j) to /(xj) + df(x.j) for each of the j = 1,...,N, with 
N an arbitrary integer and specified positions xi, x 2 , . . . ,xjy. 

Here we restrict ourselves to the study of Gaussian random fields, whose statistical properties 
are completely characterised by some power spectrum (spectral density) or its Fourier transform, 
the autocorrelation function. There are both physical and statistical arguments in favour of the 
assumption that the primordial density field in the Universe was indeed of this nature. If the very 
early Universe went through an inflationary phase, quantum fluctuations would generate small- 
amplitude curvature fluctuations. The resulting density perturbation field is generally a Gaussian 
random process with a nearly Harrison-Zel'dovich scale-invariant primordial power spectrum. But 
even while inflation did not occur, the density field /(x) will be nearly Gaussian in the rather general 
case that its Fourier components /(k) are independent and have random phases (cf. Scherrer 1992). 
The Fourier decomposition of the field at a specific location x can then be seen as the superposition 
of a large number of independent random variables that are drawn from the same distribution. By 
virtue of the central limit theorem the distribution of this field will approach normality, and (at 
least) for small N the multivariate distribution Vn is multivariate normal (Gaussian): 



exp 

V N = 



N 

[(2*)"(detM)]V* ( 2 ) 

where M _1 is the inverse of the N x N covariance matrix M, the generalisation of the variance 
a 2 in a one-dimensional normal distribution. M is completely determined by the autocorrelation 
function £(r) if the field is a Gaussian random field, 

Ma = (/(x i )/(x J )) = £(x, - Xj ) = £(1* - x,|), (3) 

Throughout this paper the brackets (...) denote an ensemble average. The last relation in equa- 
tion (3) reflects the fact that our field is a homogeneous and isotropic random process. Since we 
can consider / as an TV-dimensional column vector, we can also write the covariance matrix M in 
the convenient form 

M = (//*), (4) 

with /* the transpose of /. By taking the limit as N — ► oo with uniform spatial sampling, the 
summations appearing in equation (2) may be turned into integrals. The result 

V[f]=e- s Wv[f] , (5) 

is similar to the quantum-mechanical partition function in path integral form, where S is the action 
functional. Although there is no direct connection with quantum field theory, S will be referred to 
as the action. The expression for the action S for a Gaussian random field can be obtained from 
equation (2), 



S[f] = \ JtbuJ dx 2 /(x 1 )K(x 1 -x 2 )/(x 2 ), 
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where K is the functional inverse of the correlation function £, 



dxif(xi - x)f (x - x 2 ) = <5 D (x! - x 2 ), 



(7) 



and 5d the Dirac delta function. The measure T> [/] is most easily evaluated on a lattice, where it 
is just the product of differentials dfi divided by a normalization constant. 

Note that we use the notation V[f] to refer to an infinitesimal probability with measure V [/]; 
the probability density is exp (—£[/]). The square brackets in V[f] and S[f] indicate that these 
are functional, i.e., they map the complete function /(x) to one number. 

To compute expectation values (A) of properties (functionals) of the random field /(x), such 
as the galaxy mass distribution or the distribution of cluster shapes, one integrates the functional 
over all possible density fields /(x), weighting each by the probability from equation (5), 



This is exactly analogous to the sum over histories or paths in the Feynman path integral formu- 
lation of quantum mechanics (Feynman and Hibbs 1965). As in quantum field theory, there are 
two practical ways to evaluate cosmological path integrals: perturbation series and Monte Carlo 
integration. 

The perturbation series approach to path integrals, based on Feynmann diagrams, is limited 
to a small number of applications, as it runs into difficulties when cosmological structures become 
nonlinear. A more general way to evaluate path integrals, which is adopted here, is by Monte Carlo 
integration. By generating realisations /$ of the density field, and evaluating the corresponding 
values A[fi\ of the quantity A, the mean of these values is determined: 



The subsequent non-linear evolution is treated by performing A-body simulations of a specific 
realisation /. The central issue in this method is the need to draw samples which have a 
probability distribution proportional to exp(— S[f]) (eq. 5). 

2.2 Gaussian Random Fields: constraints 

The complicating factor in generating Gaussian random density fields subject to one or more 
constraints is that correlations couple all points of the field with all other points. Therefore, 
instead of describing the field in terms of an infinite product of one-dimensional probabilities, one 
is forced to formulate the problem using infinite-dimensional probability spaces (see eq. 5). 

The strategy followed by Bertschinger (1987) is to incorporate the set of constraints imposed on 
the density field /(x) in the action S[f], according to the definition in equation (5). A realization 
of the constrained density field is then obtained by properly sampling the resulting distribution 
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(A) = 
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function exp(— S[f]). To make clear how the constraints are incorporated in the action, we consider 
a field /(x) that is subject to a set of M constraints, 



T = {C i = C i [f;x i ]=c i ; i = l,...,M}. (10) 

The constraints are therefore imposed by forcing the field Cj[/;x], {i = 1, . . . , M), a functional of 
the field /(x) as well as a function of the point x, to have the specific value q at the position Xj. 
The constraints Cj are assumed to be linear functional. Examples of such functional are the value 
of the field itself at the point x a , the derivative of the field /(x) at the point x^, or a convolution 
over /(x) with some function g(x), 

d 

C/3[/; x /3] = ^/( X )lx 5 =C/3, (11) 



Ct[/; x 7l = y #( x 7 - x) /(x) dx = c 7 . 



The constraints C a and can be considered as particular cases of a convolution of /(x) with 
functions y a and ^ respectively, 

g a (xq - x) = (5 D (x a - x) 

5 (12) 
gpfo-x) = — 5 D ( X/3 -x). 

A broad class of constraints can be considered as such, so that a treatment of the constraints in 
the form of a convolution is not a serious restriction. In particular, we will see in section 4 that 
the expressions for the 10 constraints needed to specify the height, shape and orientation of a peak 
in the filtered density field /^(x), the 3 constraints to specify its peculiar acceleration and the 5 
constraints to specify its tidal field can all be written as convolutions over the field /(x), with the 
convolution functions g depending on the kind of constraint. 

Since we limit our fields /(x) to those that obey the set of M constraints T, the probability of 
possible realisations /(x) is the conditional probability 7'[/(x)|r], 

v[m =-7^=w\- (i3) 

The second equality follows because the constraints are linear functionals of /, so that the joint 
probability space for / and T is the same as the probability space for /. Because the constraints 
Ci are linear functionals the central limit theorem assures them to have a Gaussian probability 
distribution when applied on a Gaussian field /(x). The covariance matrix Q of the constraints' 
probability distribution can be expressed as (cf. eq. 4), 

Q = (CC t ), (14) 
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where C is the M-dimensional column vector with elements Q, and C* its transpose. The joint 
probability V [r] for the set of constraints T is therefore the following multivariate Gaussian distri- 
bution (cf. eq. 2), 



V[T] 



exp 



or, in a more concise notation, 



[(2;r) M (detQ)] 1 /2 



V[T] = exp[--C t Q- l C) V[T], 



where the measure T>[T] is defined as 



M 

n dc, <. 

i=l 



(15) 



(16) 



M 



V[T] 



[(2vr) M (detQ)] 1 / 2 



n ,/r - 



(17) 



When each field /(x) is represented by its value at A" points (e.g. in a discrete computer repre- 
sentation) we can picture the problem in a geometrical way. The fields /(x) can be considered 
as TV-dimensional vectors (/i, ■ ■ ■ ,Jn)- The constraint set T carves out an (iV — M)-dimensional 
hypersurface in this iV-dimensional vector space, consisting of all fields obeying these constraints. 
In other words, the set T is an (N — M)-dimensional hypersurface, in particular a hyperplane when 
the constraints are linear. 

The expression for the conditional probability of the field /(x) given the set of constraints T, 
V[f\T], follows after inserting equations (5), (6) and (16) into equation (13), 



V [f\Y\ = exp 



/(xi) J K'(xi - x 2 )/(x 2 ) d Xl dx 2 - CQ^C) 



V[f] 
V[V\ 



(18) 



This result shows that the constraints T are incorporated into the formalism by a change of the 
action S[f] to 



2S[f] = JI /( Xl )K( Xl - x 2 )/(x 2 ) d Xl dx 2 - C'Q-'C 



(19) 



In Appendix A it is shown that this constrained action may be written in a simple and revealing 
form, 



2S[F] = J dxi J dx 2 F(x 1 )K(x 1 -x 2 )F(x 2 ), 



(20) 



where the residual field F(x) is defined as the difference between a Gaussian field /(x) satisfying 
the constraint set V and the ensemble mean /(x) of all these fields, 



F(x) = /(x) - /(x). 
7 
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Figure 1. Illustration of the construction of a constrained random field. The field contains two peaks, an elongated one defined 
on a Gaussian scale of Ah^ 1 Mpc at [x, y] = [65.0, 65.0] h~ 1 Mpc, and a more compact one defined on a Gaussian scale of 2h~ x 
Mpc at a position of [x, y] = [35.0, 35.0] Mpc. The corresponding mean constrained field (/) is shown in the top left frame, 
to which the residual field F = / — / in the top right frame is added to obtain the constrained random field realization (/) 
shown in the bottom frames. The left frame shows the field after smoothing with a Gaussian filter with Rf = 2h~ 1 Mpc, while 
the right frame is the field after smoothing on a scale of Rf = 4/i _1 Mpc. The fluctuation field has a standard cold dark matter 
spectrum (Q = 1.0, h = 0.5). 

The conditional probability function can therefore be described as a shifted Gaussian around the 
ensemble mean field, /(x) (see Appendix A), 

7(x) = (/(x)|r)=^(x)^ 1 c„ (22) 

where summation over repeated indices is used. Thus, /(x) is the "most likely" field satisfying the 
constraints and it equals the "average density profile" obtained by BBKS. More precisely, / = / is 
a stationary point of the action: 

5 £ = for / = /. (23) 
of 

In equation (22) £j(x) is the cross-correlation between the field and the ith constraint Cj[/;Xj] 
while £ij is the (ij) th element of the constraints' correlation matrix Q, 

e i (x) = (/(x)C i ), 

(24) 

dj = (CiCj). 

If the constraints involve only the field itself at single points, like C a in equation (11), both the 
correlation matrix £jj and £j (x) reduce to the two-point correlation function £ (x) , 

e i (x) = (/(x)/(x i ))=e(|x i -x|), 
^ = (/(x i )/(x j ))=£(|x i -x j |). 

In effect, the residual field F(x) provides random noise which is added to the signal /(x), which 
is completely fixed by the imposed set of constraints T. Generating a sample /(x) obeying the 
constraints {Cj[/;xj] = c$; i = 1,...,M} therefore consists of constructing / from Cj and q 
according to equation (22), subsequently generating the noise F(x), and adding them: 

/(x) = /(x) + F(x) = -ix)£,/r, + F(x). (26) 

Notice that the residual field F is a Gaussian field because it is the difference between two Gaussian 
fields. The whole problem of constructing a constrained random field has now been reduced to a 
proper sampling of F. This is complicated by the fact that F(x) is not entirely random but subject 
to the set of M constraints Tq: 

r = {C i [/;x i ]=0; i = l,...,M}. (27) 

This follows directly from the fact that the constraints C« are linear functionals and F is the 
difference between two fields, 
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Figure 2. Linear density profiles along the central x-axis of the field shown in figure 1. On the left, the field has been smoothed 
using a Gaussian filter with Rf = 2h~ 1 Mpc. The solid line shows the constrained field (/). The dotted line is the mean field 
(/) and the dashed line the residual field F = / — /. On the right the same field, but now after filtering on a scale of 4/i _1 Mpc. 

An illustration of the sketched constrained field construction procedure, based on equation (26), is provided by figure 1. 
Note that both the original Bcrtschingcr prescription (1987) and the Hoffman-Ribak procedure (1991) are based on this equation 
(the particular realization in figure 1 has been generated with the Hoffman-Ribak code described in this paper). The fluctuation 
field in the lOOh. -1 Mpc box has a standard cold dark matter spectrum (Q = 1.0, h = 0.5) and contains two peaks of different 
shape and scale, a spherical 4<ro(2/i -1 Mpc) overdensity and an elongated 3ao(4h~ 1 Mpc) ovcrdensity. Density contour maps 
(filtered on a scale of 2h _1 Mpc) of the mean field / defined by this constraint (top left), an accompanying residual field 
realisation F (top right) and the resulting constrained field / (bottom left) are shown in slices of width 1 /20th of the boxsizc 
taken along the ^-direction. The slices pass through the centre of the box. Figure Id shows the constrained field / smoothed 
on a scale of 4/i — 1 Mpc. A good idea of the relative amplitudes of the mean, residual and constrained field in figure 1 can be 
obtained from linear density profiles through the density field. Figure 2 shows such profiles, taken along the central x-axis, 
passing through the outskirts of both peaks. The left figure corresponds to the density field at a Gaussian filtering scale of 
2h~ 1 Mpc, while the right figure has a Gaussian smoothing scale of 4h~ 1 Mpc. The dotted line is the mean field /, the dashed 
line the residual field F, and the solid line the superposition of the two, the constrained field realization /. 



Ci [f] = a [/ - /] = a [/] - a [/] = Ci - Cl = o. (28) 

This fact is independent of the numerical values {q} of the constraints V imposed on the field /(x). 
3. Sampling constrained Gaussian random fields 

Application of the construction procedure based on equation (26) requires the ability to properly 
sample exp(— S[F\) for the random field F(x). The sampling procedure forms the core of any 
constrained random field algorithm, and determines its effectiveness and reliability. The sampling is 
carried out most conveniently in Fourier space, where the action S[F] is diagonalized (appendix B), 

_ f dk |F(k)| 2 

S[F] - 7 (2^5* ' (29) 

where F(k) is the Fourier transform of the residual field F(x), 

F(x) = J T&o 5 F{k) e " kx ' (30) 

and P(k) the power spectrum of the field (see eq. 41 for the formal definition). Note that in this 
paper we adopt a different Fourier transform than Bertschinger (1987, 1992). 

In the case of an unconstrained field, for which F(x) = /(x), all harmonics -F(k) are mutu- 
ally independent and normally distributed. This makes sampling relatively easy. However, for a 
constrained field the residual field is subject to the constraints To (eq. 27), so that its Fourier com- 
ponents are no longer mutually independent. The coupling of the different Fourier modes turns the 
sampling of the action S[F] into a non-trivial problem. In earlier work we (Bertschinger 1987, Van 
de Weygaert 1991) accomplished the sampling of the residual field, carried out in discrete Fourier 
space F(kj), by means of an iterative "simulated annealing" technique. The action was sampled 
by means of a Markov chain, starting with an initial guess for the harmonics and updating them 
iteratively, each update depending only on values of the most recent estimate. After a number of 
iterations the Markov chain relaxes to a steady state with F correctly sampling the action. The al- 
gorithm used for the update is the "heat bath" algorithm, which treats the discrete set of harmonics 
F like a series of coupled oscillators in thermal contact with a heat bath of fixed temperature. The 
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heat bath generates random fluctuations in each harmonic which couple to all other harmonics. 
The fluctuations drive the system towards a state of "thermal" equilibrium in which the action is 
distributed properly. The algorithm requires 0[(M 2 + 1)N] operations to generate one independent 
realisation, where N is the number of degrees of freedom (roughly the number of grid points for 
the density) and M is the number of constraints. A disadvantage of this iterative approach is 
that as the grid density grows and as the number of constraints increases to more than a few, the 
system "anneals" so slowly that the algorithm becomes prohibitively expensive and impractical. 
Additionally, there is no unique way of deciding at which stage the system has annealed to the 
desired equilibrium. 

3.1 Hoffman-Ribak Algorithm 

The crucial observation by Hoffman &; Ribak (1991) is that the residual field F(x) has some unique 
properties which simplify the construction of a realisation of a constrained field substantially. While 
it was already known that the mean value of -F(x) is independent of the numerical values Cj of the 
constraints T, 



it had not been realized earlier that this is true for the complete probability distribution V [F\T] of 
the residual field F(x) itself (see appendix C), i.e. 



The observation that the statistical properties of the residual field F(x) are all independent of 
the numerical values c% is the key element of the Hoffman-Ribak method, rendering unnecessary a 
direct sampling from the complicated action S[F]. A particular residual field F(x) can as well have 
been sampled from the set of fields subject to the constraints T as from the fields belonging to some 
arbitrary constraint set T. The residual field -F(x) that is obtained by generating an unconstrained 
realisation /(x) of the field, and subtracting the mean field / of the constraint set Y to which it 
belongs, is therefore a correctly sampled residual field for the constraint set I\ 

These considerations lead to the following strategy for constructing a constrained realisation 
of the field /(x), consisting of five stages: 

(1) Create a random, unconstrained, realisation /(x), a homogeneous and isotropic Gaussian ran- 
dom field whose statistics are determined by the power spectrum alone. 

(2) Calculate for this particular realisation /(x) the values q of the constraints {Cj(x)| x . ,i = 
1, . . . , M}. These variables define a set of constraints, T = {c{\. 

(3) Calculate for this "random" constraint set V the corresponding mean field, using 



<F(x)|r> = </(x) - 7(x)|r> = </(x)|r> - 7(x) = o, 



(31) 



V[F\Y 1 ]=V[F\Y 2 } for all Y 1 , Y 2 . 



(32) 



/(x) = </(x)|r> = 



(33) 



(4) 



Evaluate the residual field F of the random realisation: 



F(x) = /(x) - /(x). 



(34) 
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This residual field F can also be considered the residual field of a particular realisation subject 
to the desired constraints, T. 
(5) Evaluate the desired mean field /(x), using equation (22), and add it to the residual field F(x) 
(eq. 34) to obtain a particular realisation of the desired constrained Gaussian random field 
/(x): 

/(x) = /(x)+Ci(x)^ 1 (c J --c j ) (35) 

The field /(x) constructed in this way obeys the constraints and replaces the unconstrained field 
/(x). Note that there is a one-to-one correspondence between the trial field /(x) and /(x). Fur- 
thermore, the ensemble of realisations produced by the algorithm presented here properly samples 
the subensemble of all realisations constrained by T. The algorithm is optimal because it is exact 
and involves only one realisation of an unconstrained random field and the calculation of the mean 
field under the given constraints. 

3.2 The practical implementation 

Our implementation of the Hoffman-Ribak algorithm has two important elements. Firstly, for 
reasons of convenience, all necessary calculations are carried out in Fourier space. Secondly, the 
constrained field /(x) is generated on a periodic three-dimensional lattice of side L, so that /(x) 
is evaluated on iV(oc L 3 ) gridpoints. The result can be considered to be an iV(oc L 3 ) vector 
/ = [/(x 1 ),...,/(x JV )]. 

The central equation of the Hoffman-Ribak algorithm for generating a constrained field real- 
ization /(x) is equation (35). We assume that, as in the case of the 18 peak constraints (section 4), 
the M constraints Cj[/; Xj] = Cj on the field /(x) are convolutions of the field /(x) with some kernel 
#i(x;xj), 

Ci[/;xj] = J dxi?j(x;xj)/(x) = Cj. (36) 

In the case of the peak constraints on the local density field (section 4.2) the convolution kernel is 
a Gaussian filter function or one of its first or second derivatives. 

The Fourier transforms of the field /(x) and the kernel flj(x;Xj) are defined by 

r <7k - 

'< x >= /77l3/(l<>e~* k "^ 

( £ < 37 > 



F i (x;x l )= J — -3^(k)e- kx 



(2*) d 

Consequently, Parseval's theorem yields the following Fourier expression for the constraint 
Ci[f; Xj] = Q, 



dk 

The constraint's correlation function £jj can be evaluated by using equation (38), 



/dk 
— - J fl?(k)/(k)=c i> (38) 



11 



= (Q[/;x,] C,-[/;x,-]> = (/ iT*^) /(k^ | ^ F,(k 2 ) /* (k 2 )| 



dki <ik 2 
(2^(2^0 

This immediately leads to the Fourier integral expression 



& = / 7^3^(k)^(k)P(fe), (40) 



dk 

where we have used Bertschinger's definition (1992) for the spectral density P(k), modified by a 
factor (27r) 3 owing to our different Fourier transform convention, 



(2vr) 3 P(A: 1 ) «y D (ki - k 2 ) = (/(ki)/*(k 2 )) , (41) 

with (5i)(ki — k 2 ) the Dirac delta function. Once the expression for P(k) and the Fourier transform 
Hi(k) of the constraint kernel are known, £y can be easily calculated from equation (40). 

In a similar way we obtain an expression for the cross-correlation between the field and the i th 
constraint, £j(x), 

e i (x) = (/(x)C i [/;x i ]> = ^y ^3 7( kl ) e - iki - x / ^3^(k 2 )/*(k 2 )) 

(42) 



/ 



dkl ^ (/(ki)/* (k 2 )) ^(k 2 ) e- iki x . 



(2vr) 3 (2vr) 3 

which in combination with the definition of the spectral density (eq. 41) yields the expression 

fc(x) = _/ ^3^(k)P(fc)e- k - x . (43) 

Inserting this expression into equation (35) leads to the following Fourier integral expression for 
the constrained field, 

/(x) = /(x)+^(x)^ 1 (c j -c j ) 

(44) 



dk 



(2k) 



F(k)+P(k)H t (k)^( Cj -c 



e -ik.= 



The only element left in the calculation of the constrained realization /(x) is the unconstrained 
field /(x). As was noted above, /(x) is most conveniently generated in Fourier space, where its 
Fourier components F(k) are mutually independent and Gaussian distributed. 

In practice the above expressions are evaluated on a three dimension grid of iV(oc L 3 ) gridpoints, 
and the corresponding Fourier integrals are replaced by discrete Fourier sums. Summarizing, the 
process of setting up a constrained field for a given spectrum P(k) consists of four steps. Firstly, 
the value of the constraint kernel is evaluated on N Fourier gridpoints kj, an 0(N) operation. 
Secondly, the matrix ^ is calculated by means of equation (40), with a total computational cost 
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proportional to 0(M 2 N), after which its inverse is determined, a 0(M 3 ) procedure. Subsequently, 
the N unconstrained field components F(k) are generated, from which the value of the correspond- 
ing constraint values q are evaluated using equation (38). The computational cost of the latter 
is 0{MN). Finally, the constrained field / is determined from equation (44), consisting of the 
0(M 2 N) evaluation of the products i/j(k)^ 1 (cj —cj) for all wavenumbers k, followed by a Fourier 
transform of cost 0(N log N). Thus, the total cost is 0[(M 2 + log N)N] (the cost of inverting the 
constraints is negligible because N 3> M). Although this scaling is no better than the 0[(M 2 + l)iV] 
scaling of the iterative heat bath method (Bertschinger 1987), the coefficient of proportionality is 
much smaller because no iteration is required. 

4. The peak constraints 

An important cosmological application for a constrained random field algorithm is the generation 
of an initial density field containing one or more peaks (or, equivalently, dips). A peak is identified 
as a local maximum in the density field that has been smoothed by some filter function or, more 
generally, as the immediate surroundings of this maximum. The choice of the filter will depend 
on the specific application. The scale of the peak is defined to be the characteristic scale of that 
filter function. Depending on their scale, these density peaks may be the progenitors of galaxies, 
clusters or superclusters. The constrained random field algorithm makes it possible to specify the 
height, compactness, shape and orientation of the density field in the immediate vicinity of the 
peak, while the total matter distribution can be sculpted such that the peak is subjected to a 
desired amount of net gravitational and tidal forces. In the linear clustering regime these forces 
are directly proportional to the peculiar velocity of the peak and the components of the shear at 
its location. 

Unlike the other constraints, the four quantities to describe the position and scale of the peak 
are not imposed via the algorithm described in the previous section. Rather, they are parameters 
that enter via the kernels iJj(x;Xj) (eq. 36) of each of the constraints. In addition to its scale and 
location, a peak in the smooth density field is specified by 18 constraints. The height of the peak 
needs to be specified while 3 constraints are needed to ensure that the 3 first derivatives of the 
smooth density field vanish at its summit. The 6 second-order derivatives of the density field are set 
by specifying the compactness, the axis ratios and the orientation of the peak. These 10 constraints 
together determine the density distribution in the immediate vicinity of the peak. The specification 
of the gravitational field around the peak introduces 8 additional constraints: The 3 components of 
the smoothed peculiar acceleration at the location of the peak and the 5 independent components 
of the traceless tidal field tensor. 

The constraints Cj that we use in our peak algorithm are a combination of one or more of the 
above quantities. Once a constraint has been specified an expression for the corresponding kernel 
Hi is derived (see eq. 36). In the practical implementation we derive the expression for the Fourier 
transform of Hi, iij(k). By working directly in Fourier space we save one FFT and at the same 
time guarantee a higher accuracy of the results. 

After an initial phase of linear evolution in which the Zel'dovich (1970) approximation is used, 
the further non-linear evolution of the matter distribution surrounding the peak is usually followed 



13 



by N-body simulations. It is evident that the use of the constrained random field code makes it 
possible to study the formation and evolution of these objects more systematically than possible 
with the conventional methods based on unconstrained fields. Among others, this will provide 
considerably more insight into the question of which physical parameters and processes have the 
largest influence on the fate of an object. 

In the following we will drop the explicit time dependence in our notation. The value of each 
of the quantities will be the value that the quantity has when it is linearly extrapolated towards 
the expansion factor a (with a = 1 the present epoch). The treatment in the next sections will 
be in comoving coordinates and wavevectors, while all spatial derivatives are with respect to these 
comoving coordinates. 

4.1 Peak scale and position 

Many cosmological studies have assumed that present-day nonlinear object like galaxies or clusters 
are the result of the collapse of peaks in the primordial density fields whose height exceeds some 
threshold, after having smoothed the field with a filter of a certain shape and scale. Because many 
cosmological scenarios do not posses a natural filtering scale, often an ad hoc filter has to be invoked 
to define the objects. In this paper we use a Gaussian filter because of its simplicity and smoothing 
properties. However, the formalism is equally valid for any other filter, and it is trivial to modify 
the equations (or our computer program) correspondingly. 

Although the precise relation between the Gaussian filtering scale Rg and the characteristic 
mass M p k of a particular object in the present universe is unclear — indeed, the one-to-one associ- 
ation between objects and density peaks is questioned by recent works (Katz, Quinn & Gelb 1993, 
Bertschinger Sz Jain 1994, Van de Weygaert & Babul 1994) — we can estimate a reasonable choice 
using a simple argument. The total mass enclosed by a Gaussian smoothing function with filtering 
scale Rg in a homogeneous Einstein-de Sitter universe of density p is 

M pk {R G ) = (2tt) 3/2 pR% = 4.3718 x 10 12 R% /i _1 M , (45) 

where Rg is in units of /i _1 Mpc. For example, if we take for M p k the typical mass of the core of a 
cluster, M c = 6 x 10 14 M , this yields a Gaussian filter scale of Rg ~ 4h~ 1 Mpc. Similarly, a radius 
of Rg ~ 0.6/i _1 Mpc corresponds to a mass of « 10 12 M , comparable to the mass of a galaxy with 
a luminosity equal to L* if f2 = 1. 

The use of the filter function Wg serves a twofold purpose in our peak constraint algorithm. 
In addition to defining the scale of the peaks in the density field p(x) it is also of vital importance 
in the derivation of the kernels Hi of each of the constraints. The expressions for these kernels are 
found by using the fact that the peaks are maxima in the filtered density field /g(x), 

/ g (x) = J dyf(y)W G (y,x), (46) 
where /(x) is the density contrast field, 
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/(x) = (47) 

(in this equation p is the average density of the Universe). This convolution integral is equivalent 
to the Fourier integral 

feW= / — 1I /(k)H"(k;x) 

/■ 4 - < 48 » 
= /(W /<k)W '* (k>e ^ 

where IU(k;x) and VT(k) are the Fourier transforms of VFc(y,x) and Wg(x, 0). In the case of a 
Gaussian filter, 

= (2^72 (- i ^ L ) ' (49) 

#"(k;x) and W(k) are 

W^(k) = e- k2R a/ 2 and W(k, x) = W(k) e ik ' x = e^c/ 2 e ik * . (50) 

Note that the position x of an object causes a phase shift k • x with respect to an object that is 
situated at the origin, 0. 

4.2 The local density field 

Locally, the density field around a peak at position x^ can be described by the second order Taylor 
expansion of the density profile /g( x ) around the peak, 

/g(x) = /o(xd) + - ^ q (x d ) (xi - x dti )(xj - x dJ ) . (51) 

i,j=l 1 

In this expansion we have used the fact that the three first derivatives of the field /g(x) at the 
location of the local maximum, x<f , are equal to zero. The equation shows that the requirement that 
the smoothed density field /g( x ) has a maximum of a certain height, shape, orientation at location 
Xrf translates into constraints on the value of the smooth density field /g at x^, on its gradient 
V/g and on the second derivative tensor of the field, VjVj/c- This implies that 10 constraints are 
required to fully specify the local density field around a peak. Also note that the quadratic part 
of equation (51) should be negative definite if /g(x^) is a maximum. Consequently, the isodensity 
surfaces fc = F around the peak are triaxial ellipsoids, whose orientation and size depends on the 
value of the second derivatives of fa- 

The first constraint is the height of the peak, /g(x^). Usually it is expressed in units of the 
variance ao(Rc) = (/g/g) 1 ^ 2 of the smoothed density field, 

/ G (xd) = u c a (R G ) , (52) 
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which in combination with the convolution expression for Jq in equation (46) yields the following 
expression, 



r dk ~ 

J f&o 1 /(k) w * {k) e " kXd = Uc(To{Rg) ■ (53) 

Consequently, the corresponding constraint kernel Hi(k) (see eq. 38) and constraint value c\ arc 
given by 

&i(k) = W(k) e ik ' Xd , d = iy c a (Ra) . (54) 

For reasons of clarity and convenience a compilation of the kernels of all peak constraints is given 
in appendix F. 

Three additional constraints are obtained from the extremum demand that the first order 
derivatives of fc should be at the peak position x^, 

^(xd) = 0, i = 1,2,3. (55) 

The Fourier expression for the gradient V/c(xd) is obtained by partial differentiation of the inte- 
grand in the convolution equation (48), 

This yields the following constraint expressions, 



/ 



dk ,k 



ikf(k)W*(k)e~ lk ^ = 0. (57) 



The corresponding kernels ^(k), ^(k) and i?4(k), and the constraint values C2, C3 and C4 are 
therefore 

Hj(k) = ikW(k) e ik Xd , Cj = , (58) 

where j = 2, . . . , 4 and the corresponding I = j — 1 (also see appendix F). 

Finally, there are six constraints that correspond to the shape, compactness, and orientation 
of the density field around the peak. Because the density field in its vicinity is ellipsoidal (see 
appendix E), its shape is fully characterized by the two axis ratios a\i = a\jai and 013 = 01/03. 
The quantity that describes the compactness, or steepness, of the density profile around a peak is 
the Laplacian V 2 /g(x^). Usually this Laplacian is expressed in units of U2{Rg) = (V 2 /gV 2 /g)^ 2 
(see appendix E), 



V 2 /c(xd) = -x d (T 2 {R G ) 
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(59) 



The minus sign in this definition of Xd is introduced in order for Xd to be negative in the case of a 
dip and positive for a peak. The orientation of the peak with respect to the coordinate axes is most 
conveniently specified by the three Euler angles a, (3 and ip. The corresponding transformation 
matrix Aij is given by, 



cos a cos ip — cos (3 sin a sin ip sin a cos ip + cos (3 cos a sin ip sin (3 sin ip 
A = | — cos a sin ip — cos (3 sin a cos ip — sin a sin ip + cos /? cos a cos V> — sin (3 cos ) • (60) 
sin (3 sin a — sin /3 cos a cos (3 

The above six quantities (a±2, ai3, x^, a, f3 and ■0) constrain the six second order derivatives of fa 
via the combination (see appendix E for a derivation), 

~ a = ~ X! ^kAkiAkj , i,j = 1,2,3, (61) 

where ^4^- are the elements of the orientation matrix (eq. 60), and the Aj are the eigenvalues of the 
matrix — VjVj/c. The values of A, are obtained from the axis ratios a±2 and 013 of the isodensity 
ellipsoids around the peak, as well as from the steepness of the density profile, Xd, via the relations 

x d a 2 (RG) . , 2 \ \ 2 /« x 

^1 — 77— 2 2\ ' A 2 — Aia 12 , A3 — Aia 13 , (b^J 
i 1 + a 12 + a 13i 

A Fourier expression for the second order derivatives of /g( x ) is obtained by double partial differ- 
entiation of the integrand of the convolution integral (eq. 48), 



/ 



(27T) 

so that we find 



k=i 



H t (k) = -kikjW(k) ^ , Q = - E X kA ki A kj , (63) 

k=l 

for the kernels Hi(k) and constraint values q, with I = 5, . . . , 10 and i,j = 1, . . . , 3 (see appendix F 
for the correct numbering). 

4.3 The local gravitational field 

The peak constraints that were introduced and discussed in section 4.2 describe the density field 
in the immediate surroundings of the peak. Of more fundamental importance to the dynamics of 
a region of space are constraints on the gravitational potential perturbations. The local poten- 
tial perturbation 0(x) is the weighted sum of all density perturbations throughout the universe. 
Constraints on the local potential therefore have immediate repercussions for the global matter 
distribution. Since we wish to neglect the potential fluctuations on scales smaller than the objects 
we are interested in, we consider the smoothed potential perturbation field 4>g, 



g (x) = J dy <f>(y) W G (y, x) . 
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(64) 



In our implementation we use a Gaussian function for the filter Wc;(y,x), as in the case of the 
density field. 

It is physically appealing to impose constraints on the potential (f> via constraints on its deriva- 
tives, in particular the gravitational acceleration and the tidal field. The peculiar gravitational 
acceleration g at the position x is 

1 d(av) 1 , , 
a dt a 

where a is the cosmological expansion factor and v(x, t) the peculiar velocity of the patch of matter 
at physical position r(i) = ax(t), 

v=%-Hr. (66) 

A first-order Taylor expansion of the gravitational field around the peak shows that the dynamical 
state of the patch of matter in its immediate neighbourhood, on scales larger than the filter scale 
Rg, is completely specified by the bulk acceleration gc(x^) = — V^c/a, the divergence V • gc and 
by the traceless (comoving) tidal tensor Ea,ij, 



9G,i( x ) = 9G,i( x d) + a Yl {^a~( V ' SG)(xd) Sij - (xj - x dJ ) , 



(67) 



where Sy is the Kronecker delta and Eo,ij is the trace- free part of —dgG,i/drj = d 2 (f>G/dridrj (note 
that here we choose to use physical coordinates rj, since we are dealing with physical quantities), 

j? 1 I dg G ,i . dg GJ \ 1 if d 2 4> G 1 r7 2, x \ / ^!Q ^ 

EG ^ = ~Ya {^7 + -dxjj + 3^ (V • gG) = 7 \dx~dx- ~ 3 V • (68) 

The divergence V • gc/o is the component of the gravitational field corresponding to pure radial 
infall into (or outflow from) the peak. Through the Poisson equation this quantity is directly 
proportional to the local density perturbation /g( x )> 

^ = -\^c = -l^H 2 / G (x) . (69) 



The expression for the constraint on V • gc/a is therefore equivalent to equation (53), except for 
the proportionality constant 2>QH 2 /2 in both constraint kernel Hj and value Cj. From the above 
equation we can also easily infer the relation between the Fourier components 0c(k), 

<fc(x) = / -^3^(k)e- lk ' x , (70) 



„ (2vr) 3 

and the Fourier components /(k) of the density field, 



fc(k) = -^vi/(k)r(k). (7i) 
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The first 3 constraints on the gravitational field therefore concern the peculiar gravitational accel- 
eration at the position of the peak itself, gG(x<f). It is useful to specify it in units of the dispersion 
of the gravitational acceleration of peaks, cr g>p k(RG) = (gG,pfc ■ EG,pk), 

9G,l(*d) = 9l<Tg,pk{RG) , 1 = 1,..., 3. (72) 

The dispersion of the peak accelerations is less than the overall dispersion a g of the acceleration in 
the field. This lowering of the acceleration of peaks compared with that of field points is caused by 
the extra acceleration associated with the infall of field points onto the peaks. We can infer that 
(see section 4.4, eqns. 101 and 106) 

<7g,pk = °g = Vg \A — lv > ( 73 ) 

where a g (Rc) and j v are given by 

3 a 2 
a g (R G ) = -QH 2 g-^Rg) and lv = — 2- , (74) 
Z cr_lC7l 

with Uj{Ra) the spectral moments, 

«](Rg) = J ^5 P(k) W(k) A*" . (75) 

The Fourier expressions for the 3 components of the bulk peculiar acceleration gci^-d) can be 
derived from equation (65) and (70), 

Inserting equation (71) and (72) leads to the constraint equations, 

From this we find the corresponding constraint kernels Hn(k), H\2(k) and #13 (k), and the con- 
straint values en, c\i and C13, 

Hj(k) = |f W(k) , a = ~ 9 i a _ l{ R G ), ( 78 ) 

with j = 11, . . . , 13 and I = 1, ... ,3 (see appendix F). Evidently, instead of specifying the constraint 
values cj as gi it is also possible to do this directly in the appropriate physical units (e.g. km/s 2 ). 

Five additional constraints are needed to characterize the tidal field around the peak. This field 
is described by the traceless (comoving) tidal tensor Eq^ (eq. 68). Within an arbitrary system 
of reference this tidal tensor is most conveniently expressed in terms of its eigenvalues and units 
vectors, 
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() ' 0(; -\y 2 <P G 5 lj ) = Y j £ k T kl T k3 , *,j = 1,2,3. (79) 




dxidxj 3 

The elements of the matrix T k i are the components of the various eigenvectors of the tidal tensor, 
whose directions are characterized by the 3 Euler angles a^, 0e, and ipE {Tki are given by equa- 
tion (60), with OLE-, Pe and ipE replacing a, (3 and ip). In an initial random density field there is 
a strong correlation between the tidal tensor and the mass tensor Qj = VjVj/ (see section 4.4, 
eq. 106). In the case of peaks this translates into a strong tendency of the tidal tensor to align 
itself along the principal axes of the mass ellipsoid. In the specification of the initial tidal field 
it is therefore often useful, and physically sensible, to express its elements with respect to the 
reference system defined by these axes. We denote the corresponding transformation matrix by 
Tku which is defined through equation (60) by the 3 corresponding Euler angles q.e, Pe and ^>e- 
If the orientation of the peak itself with respect to an arbitrary reference system is specified by 
the transformation matrix Aki(a, (3,ip) (see eq. 60), then the tidal field's transformation matrix T 
within the same system is the matrix product of T with A, 

3 

Tki = J]] Tk m A m i , (80) 

m=l 

The magnitude of the tidal field in the directions of the principal axes of the tidal tensor is 
given by the eigenvalues Si, £2 and £3. Because E G ^j is traceless, i.e. J^£k = 0, it is sufficient 
to specify two eigenvalues. To get an idea of the right order of magnitude it is usually useful to 
specify £k in units of <je, the dispersion of the off-diagonal elements of the tidal tensor E G ^j (see 
section 4.4, eq. 101), 



3 1 — 2 2 

cte(Rg) = ^H 2 a (R G ) \ —L W ith 7=—, (81) 

so that £k = £k cte(Rg)- An elegant and convenient parameterization of the diagonalized Eq^j in 
terms of two quantities e and w was introduced by Bertschinger &; Jain (1994), 

E G ,ij = diag [£ 1 ,S 2 ,S 3 ] = QH 2 e (1 + f G ) Q^w) , (82) 
with the one-parameter traceless matrix Qij defined by 



Qijiw) = diag [Qi, Q 2 , S3] = diag 

This matrix representation turns out to be useful when considering the Lagrangian equations of 
motion of a patch of matter (Bertschinger & Jain 1994). It is particularly convenient because all 
possible eigenvalues of E G ,ij are obtained by qQij(a), with q G [0,cx>) determining the magnitude 
of the tidal field, and a G [0, ir] the relative strength of the tidal field along the three principal axes. 
The 5 constraints, for Ec,n, E g ,22, ^g,12, -E-g,13 and E G ^z, therefore have the form 



W + 27T 

cos I I , COS 



w — 2tt 



( w 

COS — 

V 3 



(83) 
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3 

E G ,ij = ea E (R G ) ]T Q k (w)T ki T kj , (84) 
fc=i 

with = (1, 1), (2, 2), (1, 2), (1, 3) and (2, 3). Note that here we have expressed e in units of <te, 
i.e. £ = e Q(w). In addition, we have assumed that the fluctuations are linear, so that the factor 
ef G can be neglected. For the generation of initial conditions, our primary interest, this assumption 
is not a serious restriction. 

The Fourier components E G jj()c) of the tidal tensor, 

E °' 3 = J (2^^' (k)e " k ' Xd ' (85) 

can be easily found from the definition of E G ^ in equation (68) and subsequent differentiation and 
insertion of equation (71), 

EqM = {^p- - W*(k) /(k) . (86) 

This leads to the following tidal field constraint expressions: 

/ {I qr2 " k 6 *) ^* (k) e ~ ikxd } m = i<te{rg) P ' ^ 

which yields the corresponding 5 constraint kernels and values, 



fl,(k) = - W(k) , q = e hlH 2 a (R G ) ^ ]T Q k (w)T ki T kj , 

(88) 

with I = 13, ...,18 and = (1, 1), (2, 2), (1, 2), (1, 3) and (2,3). Alternatively, instead of ex- 

pressing the tidal constraints via the two quantities e and w (or £) and the 3 Euler angles oe, 
Pe and ipEi we can evidently specify the values for E G ,ij directly, either in corresponding physical 
units or in units of cfe(R G ), 

3 ll — 2 

E G ^ = hj ^H 2 a (R G ) y . (89) 

In the linear regime an analogous, and for some more familiar, way of describing the dynamics 
of a patch of matter is in terms of the peculiar velocity. This is possible because for growing mode 
linear perturbations the peculiar velocity v is directly proportional to the peculiar gravitational 
acceleration g (Peebles 1980), 



(90) 



where JF(O) Q, - 6 . It is convenient to write the smoothed peculiar velocity vg around the position 
x<i of the peak in terms of the bulk motion vg(x^), the divergence V • vg/a, the shear aij and the 
vorticity u>ij, 

«G,i(x) = v G ,i(xd) + • v G )(x d )<% + <Tij(xd) +o;y(x d )| (x^ - x dJ ) . (91) 

The shear is the trace- free symmetric part of dvc,i/drj, 

1 f 9w Gii <9v Gii 1 1 
^ = ^{^- + ^f)-^ (V - VG) ^' (92) 

while the vorticity Uij is the antisymmetric part, 



J_ [ dvc,i _ dvc,j 
2a 1 3xj 



(93) 



Because does not have a gravitational origin, it is an irrelevant quantity as far as constraints on 
the density perturbation field are concerned. Moroever, it can be shown to remain zero whenever 
there is no primordial vorticity (Peebles 1980). From this we can infer that constraints on the 
peak velocity v G (x^) are therefore equivalent to equation (77), except that the factor 3QH 2 a/2 in 
both the constraint kernel Hj and constraint value Cj has to be changed into HaJ 7 ^). Also, the 
constraint on the divergence V • v G /a is equivalent to the constraint on /g( x ) (eq. 53), except for 
a factor HJ 7 ^). Finally, we see that the relation between the shear ac,ij and Ec,ij is 

E = AnH 2 ^— . (94) 

We should, however, not fail to appreciate that in the nonlinear regime the simple relation between 
v and g breaks down. In that case the same basic physical relationships between the acceleration 
g, the potential <j) and the density p remain valid. As this is not true for the velocity v, it is 
fundamentally preferrable to impose constraints on the gravitational field instead of the velocity 
field. 

4.4 Probability of peak constraints 

In principle, in the case of a Gaussian field any set of numerical values for the 18 peak constraints 
per peak is possible, regardless of how small the probability of the occurrence of such peaks would 
be. This is a consequence of the ability of the Hoffman-Ribak constrained random field method to 
generate realizations for any arbitrary set of values for the imposed constraints. In order to prevent 
the generation of unlikely circumstances it is therefore necessary to control or have an estimate of 
the likelihood of the constraints. The corresponding probability distribution of the constraints is 
given by equation (15). A good measure of the likelihood can be obtained by calculating the x 2 -, 

M 

X 2 = £ CiiQ-^jCj = a&cj. (95) 
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The probability that for this constraint set x 2 has this value or higher can then be directly calculated 
from rg(M/2, x 2 /2), where Tq is the incomplete gamma function. As a rule of a thumb, the 
constraint set can be considered to represent manifest unlikely conditions, if the x 2 P er degree 
of freedom, x 2 = X 2 '/M , differs significantly from unity. Note that the computational cost of 
evaluating x 2 1S negligible as the inverse of the constraint-constraint correlation matrix ^ = (CjCj) 
has already been calculated as part of the construction procedure (eq. 40 and 44). 

A full expression for x 2 or > even better, the full probability distribution in terms of the 18 
constraint quantities can be obtained by evaluating the expression in equation (15), following the 
treatment presented in Appendix A of BBKS. Following the discussion in the previous sections the 
density and gravity field in and around an arbitrary point x in a Gaussian random density field 
/(x) can be characterized by 18 parameters 



T = {",Vi,V2,V3,(i, (2,(3,(4,(5, (a, 9i, 92, 93, E 1 ,E 2 ,E 4 ,E 5 ,E 6 } , (96) 

with / = voq the value of the field at x, Vj/ = rji the first derivatives of the field, and £4 
the six independent components of the tensor Qj = VjVj/ (where A = 1,2,3,4,5,6 refer to 
the ij = 11,22,33,12,13,23 components of the tensor). In addition, g« = — Vj^> is the peculiar 
gravitational acceleration, while Ea are the five independent components of the traceless tidal 
tensor E^ = VjVj> - jV 2 0<% (with A = 1,2,3,4,5,6 referring to the ij = 11,22,33,12,13,23 
components of Eij). 

The probability V(T) that at the position x the field has the specified values for these 18 
quantities is specified by a joint Gaussian probability distribution, for which a reasonably insightful 
expression can be found by reducing the corresponding 18 x 18 covariance matrix Q = {yiUj) into 
a block diagonal matrix of 9 2 x 2 blocks. This is achieved by transformation of the set of variables 
{d,(2,(3,Ei,E 2 } into a new set {x,y, z, E y , E z }, 

Ci + C2 + C 3 C1-C3 Ci - 2C2 + c 3 

x = , y = , z = , 

02 2d2 2(72 

(97) 

E\ — E3 Ei — 2E2 + E3 

Ey ~ g , E z - . 

V(T) is then given by 



V(T) = Ae~ Q/2 dvd 3 7]dxdydzd(4d( 5 d(Qd 3 gdE y dE z dE 4 dE 5 dE 6 , (98) 

with 



A 



3 6 5 5 



1024 vr 9 (1 - 7 2 ) 3 (1 - 7 2 ) 3 / 2 (^tf 2 ) 4 a\ a of 4 



(99) 



and 
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18 

2 (x — x*) 2 2 . r 2 i 3v ' V . I^Ca , nnr,\ 

= v H 2 h 15y + 5z H ^ h > — — + (lUUj 

3(ff~g.) 2 , (gy-^) 2 , {E-Etf *{EA-Etf 
where <7 g and ffg are defined by 



o g ee (^i?V_! yi - ^, o E ee fy.H^a^L-t. (101) 

while the various coupling quantities x*, g*, E*, E* and £7^ are defined by 



x* = 7 v , 

Z <7l 



3 3 (102) 

* Jy = 1 y(-SlH 2 )a {i , E* z = 7 z(-nH 2 )a , 

E* A = 7 £nH 2 )^( A , A = 4,5,6. 

2 (72 

(for the definitions of 7, j v and the various <jj see eqns. 81, 74 and 75). In the case of a peak we 
can further reduce the expression for Q. Evidently, fj = 0. In addition, we can use the fact that 
Q should be independent of the orientation of the mass ellipsoid around the peak, expressed by its 
Euler angles a, f3 and tp. We can therefore discard the orientation term l^Cl/ "! an d redefine 
x, y and z in terms of the eigenvalues A« of (— Qj), whose relation to the axis ratios of the ellipsoid 
are given in equation (62), 

_ Ai + A2 + A3 _ Ai — A3 _ Ai — 2A2 + A3 , , 

X ~ (J 2 ' V ~ 2a 2 ' Z ~ 2a 2 ' [ ° 6) 
Furthermore, we restrict the ordering of the eigenvalues to Ai > A2 > A3 > 0. The condition that 
A3 > is equivalent to the demand that Qj has to be negative definite, which, together with the 
constraints Vj/, is necessary and sufficient to have a local density maximum, a peak. We then find 
for the complete probability "P(T) that at an arbitrary position x there is a peak with a height 
/ = vao, a shape characterized by the parameters x, y and z, an orientation specified by the Euler 
angles a, [3 and i/j, an acceleration g, and a tidal field described by the parameters E y , E z , E4, £5 
and Eq, or, rather, that there is a peak with these parameters in the specific infinitesimal ranges 
around these values, 

V(u,x,y,z,a,P,'>p,g,E y ,E z ,E i ,E 5 ,E 6 ) 

= A\y(y 2 - z 2 )\ sin (3 e~® /2 dv d 3 r/ dxdydz dadfidip dg dE y dE z dE 4 dE 5 dE 6 . 

In this equation the constant A is 
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A = 



3 7 5 5 



256 vr 9 (1 - 7 2 ) 3 (1 - 7 2 ) 3 / 2 {UlH^f a a\ 



(105) 



and 



1-7 <r ff o-B 3cj £ A=4 



where -Ey, E4, E5 and £6 are specified with respect to the principal axis of the mass ellipsoid 
(see section 4.3, eq. 80 for the appropriate transformations). In particular, this implies that E* A = 
for A = 4, 5,6 (eq. 102). Also note that equation (100) is therefore an expression of the fact that 
in an initial random density field the tidal field has a strong preference to align itself along the 
principal axes of the the mass tensor Qj. In particular, for a peak this implies that the strongest 
tidal force tends to be directed along its smallest axis (the one with the highest eigenvalue \). 
As it explicitly takes into account this strong correlation between the initial mass quadrupole and 
the tidal field at the position of the peak, the reference system defined by the mass ellipsoid is 
therefore the most natural one to specify the initial tidal forces. The resulting expression for Q in 
the above equation is essentially the one for the x 2 °f the imposed constraints, once it is scaled 
to the appropriate filter and filter scale Rq by means of j(Rg), ^ v (Rg) and the various spectral 
moments aj(Rc)- 

The probability distribution in equation (104) is the one for having, at some arbitrary field 
position, a peak with the required physical parameters. Often, however, we are more interested 
in the more specific question what the probability V p k is that a peak at an arbitrary position has 
these imposed constrained properties. To evaluate this we need to determine the (comoving) number 
density of peaks with the constrained parameters, which can be done following the prescription in 
BBKS. To obtain V p k this specific number density has to be divided by the total comoving number 
density of peaks, n p k, whose value is (see BBKS), 

Hpk = J^JfRl = °-° 16 ^ 3 ' With R ^^V 2 - (107) 

From this we can derive the probability that a peak has a height uao, shape parameters x, y and 
z, an acceleration g and tidal tensor components E y , E y , E4, E^ and Eq. Since the orientation of 
the peak is here a less relevant quantity we integrate over the Euler angles a, (3 and ip. This can 
be done without further complications since Q is independent of these Euler angles. Note that this 
automatically implies that the tidal tensor components are specified with respect to the principal 
axes of the mass ellipsoid. We then obtain the following expression for V p k, 

P p k{v, x, y, z, g, E y , E Z ,E4, E^,Eq) dv dxdydz dg dE y dE z dE^dE^dE^ 

-n/o (108) 
= B&(x,y,z) F(x,y,z)e w/ dv dxdydz dg dE y dE z dE 4 dE 5 dE 6 . 

In this expression Q is given by equation (106), while 
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F(x,y,z) = y (yi-z 2 ){x-2z)[(x + z) 2 -(?>y) 2 ]. 



(109) 



In addition, the function Q(x,y,z) is defined such that its value is 1 when the peak constraints in 
the (x, y, z) domain are satisfied, and otherwise. These constraints are that y > z > —y, y > 0, 
to obtain the correct ordering of the eigenvalues Aj of Qj , and (x + z — 3y) > so that the smallest 
eigenvalue A3 is positive and we indeed have a peak. The constant B is given by 



This can be easily extended to other conditional peak probabilities, e.g. the chance that a peak 
of height vctq has the required parameters. However, calculating the involved expressions quickly 
becomes a very elaborate procedure. 

In the above we have mainly concentrated on the probability of constraints imposed on one 
particular peak, with the intention of providing insight into how the different constraints interrelate 
and to get an idea of the expected order of magnitude of each of the constraints. However, as we have 
seen earlier, our code allows to provide constraints on many different peaks, at different positions 
and scales. Giving analytical expressions for such constraints would quickly become a cumbersome 
and elaborate affair, due to the introduction of spatial correlations in the random field. However, 
via equation (95) the numerical value of \ 2 for these constraints can be easily computed, providing 
a good idea of their likelihood. 

5. Realizations and Applications 

The formalism developed in the previous sections allows the generation of a large variety of initial 
conditions. In this section we will visualize the procedure by providing some practical examples. 
The versatile and non-local nature of the formalism has already been emphasized in figure 1 (sec- 
tion 2), illustrating the construction of a density field that is constrained to have two peaks of a 
different scale, a different shape, and at different positions. Although for the construction of the 
density field around one central peak a few other equally or even more efficient methods have been 
developed (Binney & Quinn 1991, Bond & Meyers 1993), mainly based on a multipole expansion 
of the field, their efficiency breaks down if the constraints are imposed at more than one position. 

5.1 Field-constraint correlations 

At the core of our construction procedure is the superposition of a mean field / and a properly 
sampled residual field F (see eq. 26). The mean field / is effectively the superposition of the 
field-constraint correlation fields £i(x) (eq. 22), each weighted by a factor J2j £ij lc j- 

Figure 3 shows the mean field and five of the composite field-constraint correlation fields for the 
set of peak constraints described below. These realizations are generated in a periodic 100/i _1 Mpc 
box. As for figure 1, the power spectrum of the random field fluctuations is the standard cold dark 
matter spectrum of Davis et al. (1985) with ^CDM = anc ^ h = 0.5 — normalized such that 
a(8h^ 1 Mpc) = 1.0 at a = 1, the present epoch (Davis & Peebles 1983). The panels in figure 3 



B = 
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(110) 



16 vr 5 (29 - 6^6) (1 - 7 2 ) 3 (1 - llf /2 (\vH 2 f a\ a 



26 



Figure 3. The mean field / (top left panel) and five of the composite field-constraint correlation fields £i(r) of a set of constraints 
(see text). The panels contain the contourmaps of the density (top left, contour spacing 0.5) and correlation values (other 5 
panels) in a 5h^ 1 Mpc slice through the center of a 100fe — 1 Mpc box. The spectrum of the field is the standard CDM spectrum. 
The 5 field-constraint correlation functions are a) top middle: (f fo) (contour spacing 0.1), b) top right: (/V x /g) (contour 
spacing 0.03) c) bottom left: (/ V£/ G ) (contour spacing 0.02), d) bottom middle: (/ vg,x) (contour spacing 0.03) and c) (/ a xx ) 
(contour spacing 0.015), with fa the value of the smooth density field, vq,x the x-component of the peculiar velocity and a xx 
the ra-component of the shear tensor, all evaluated at the center of the box. 

contain contourmaps of the density and correlation values in a 5/i _1 Mpc slice centered halfway in 
the simulation box, each of the maps being smoothed on a Gaussian scale of 4/i _1 Mpc. 

All the constraints are defined on a Gaussian scale of 4/i _1 Mpc. A triaxial peak, with axis 
ratio 10 : 9 : 7, of height fa = 3oo and local density field curvature V 2 fc = (x)(T2 ~ 3.481o"2 
is positioned at the center of the box. Its major axes are slightly oriented with respect to the 
coordinate axes of the box. In addition to these local constraints, there are constraints on the local 
gravity and tidal field. Because we limit ourselves to growing mode linear perturbations we specify 
these constraints in terms of the peculiar velocity and shear. The total peculiar velocity of the peak 
is 1145 km/s, towards a direction 26.6° "north" of the positive x-axis and 22.6° out of the x — y 
plane, in the positive z-direction (note that the specified numerical values of the constraints are 
the linear extrapolations to the present epoch, a = 1). This corresponds to a value of 2.00 times 
the velocity dispersion of peaks on a scale of 4/i _1 Mpc, or 1.66 times the velocity dispersion of an 
average field point on this scale (the lower value of peak velocities in comparison with the velocity 
of field points is due to the extra component corresponding to accretion onto the peaks). The shear 
tensor at the location of the peak is orientated so that the off-diagonal terms are zero. The diagonal 
term in the x-direction, a xx , has the largest magnitude and is positive (dilation), while a yy and a zz 
are equal and negative (contraction). For illustrative purposes we have chosen a rather extreme 
value for the magnitude of the largest element of the shear tensor: lOOkm/s/Mpc on the scale of 
4.0/i _1 Mpc, ps 6.8 times the dispersion (« 14.5 km/s/Mpc) for the diagonal shear components for 
peaks (Bond 1987). A good idea of the order of magnitude of these shear tensor values is obtained 
by comparison with the value of the expansion scalar for the 3o"o peak, V • vg = — 142.47km/s. 

The specified constraints can be easily recognized in the resulting mean field, in the top left 
panel of figure 3 (contour spacing 0.5, and the positive value solid contours separated from the 
negative value dotted contours by the thick solid line corresponding to / = 0). The contours around 
the center of the box clearly reveal the presence of the elongated peak, oriented with respect to the 
coordinate axes. The global density field in the box reflects the gravity and tidal field constraints. 
The source of the motion of the peak is the concentration of mass in the upper righthand quarter of 
the frame, while the clearly discernable quadrupolar component in the matter distribution induces 
the tidal field. To understand how the different constraints conspire to produce this mean field 
it is quite revealing to study the individual field-constraint correlation functions £j(x). The five 
illustrated correlation functions £«(x) are the correlation of the field /(x) with (1) the value of 
the smoothed density at the peak position Xj, /c(xj) (top middle panel), (2) the value of the first 
derivative of the smoothed density field V x /g(xj) (top right panel), (3) the value of the second 
derivative of the smoothed density field V^/c(xj) (bottom left panel), (4) the peculiar velocity 
VG,x{^-i) (bottom middle panel) and (5) the shear component <J xx (~x.i) (bottom right panel). 
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Figure 4. The variance of constrained random field realizations. The mean field / and four different field realizations of a set 
of constraints (see figure 3) arc shown in the top left panel and the middle and right rows respectively. The panels contain the 
density contourmaps in the 5h~ 1 Mpc thick central slice in a 100/i~ 1 Mpc box. The contour spacing is 0.5. The bottom left 
panel is the contourmap of the value of the variance of the field realizations inside the slice, running from 0.0 at the centre to 
(j q ~ 0.95 at the edge of the box (contour spacing 0.05). 

The first correlation function (top middle panel) is spherically symmetric, with a value of 1.0 
near the centre, and radially decreasing to a value of 0.0 at the outer contour (contour spacing is 
0.1). Effectively, this correlation function is the convolution between the field correlation function 
£(x) = (//) and the Gaussian filter function defining the scale of the constrained object. In a 
similar fashion we can consider the second correlation function (top right panel, contour spacing 
0.03) to be the convolution of the correlation function £(x) with the first derivative of the filter 
function. This introduces the anisotropy along the x-axis, with, within distances comparable to 
the correlation radius, negative values on the left side of the peak and positive values to the right. 
Further outward the correlation function £(x) becomes negative, resulting in the sign reversal of 
£i(x). The third correlation function (bottom left panel, contour spacing 0.02) is essentially the 
convolution of the field correlation function ^(x) with the second derivative of the Gaussian function, 
d 2 Wc/d 2 x. Because this derivative has two zero-points along the x-axis we see a negative value 
near the centre, changing to positive on both sides of the centre. 

The correlation functions corresponding to the velocity and shear constraints display familiar 
patterns. The function corresponding to the constraint on the peculiar velocity in the x-direction, 
v x , (middle bottom panel, contour spacing 0.03) is a dipolar function centred on the position of 
the peak, with positive values to the righthand side of the peak (in the x-direction) and negative 
values to the left. This is evidently related to the fact that such a dipolar matter distribution would 
produce a net gravitational acceleration, and corresponding peculiar velocity, in the x-direction. 
In addition, we see that the constraint on the shear component a xx results in a clear quadrupolar 
pattern of the correlation function (right bottom panel, contour spacing 0.015). As in the case of 
the velocity-field correlation function, this is related to the non-zero tidal force xx-component that 
would be produced by such a quadrupolar mass distribution. 

Superposition of the complete set of these correlation fields £i(x), with the appropriate weight 
factors, proportional to the corresponding constraint values Cj, produces the mean field in the top 
left panel. Comparison of the different panels in figure 4 shows that several of the correlation 
function patterns can indeed be recognized in the mean field. 

5.2 Variance of realizations 

After having constructed the mean field / (former subsection) , we need to add a properly sampled 
residual field realization (eq. 26). Figure 4 provides an idea of the possible variations between the 
residual field realizations and, specifically, the resulting full field realizations. In addition to the 
mean field illustrated in the the top left panel, four different realizations are shown in the middle 
and right row of panels. All these panels are density contour maps (contour spacing 0.5) in the 
same central 5/i _1 Mpc thick slice used in figure 3. From the four field realizations we can infer 
that, for example, the mass concentration to the right, responsible for the peculiar motion of the 
peak, can vary substantially in position, shape, size and substructure. Moreover, the morphology 
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and distribution of mass clumps inside the band of matter along the x-axis, main contributor to 
the specified shear, displays an even larger variation, in particular at large distances from the peak. 

An analytic expression for the variance of the residual field at any position x follows immediately 
from the independence of the residual field distribution function from the numerical values of the 
imposed constraints (eq. 32), see Appendix D, 

(F 2 (x)|r) = al - -ix^'-ix). (Ill) 

with 

4 = (/ 2 (x)), (112) 

the variance of the density field (recall that both / and F have zero mean). The expression in 
equation (111) shows that (F 2 (x)|T) is dependent on x, and therefore implies the residual field 
(F 2 (x)|T) is neither homogeneous nor isotropic. Note that because F{x) is a Gaussian random 
field its distribution functional V [F\T] is completely specified by the variance (F 2 (x)|T). 

The lower left panel shows a contour map of the variance field corresponding to the constraint 
set in the example. Notice the perfect spherical character of this variance field, increasing radially 
outward from the position of the peak, where it is equal to 0.0, to the general field value <jq ~ 0.95 
(contour spacing 0.05). At first sight this might seem counterintuitive, as most of the applied 
constraints are non-isotropic. However, from equation (111) we see that (F 2 (x)|r) involves a 
product of all field-constraint correlation functions £i(x), independent of the actual numerical values 
Cj of the constraints. In our example all 18 peak constraints have been specified. This means that 
the anisotropy introduced by e.g. the dipole distribution corresponding to the v x constraint gets 
fully compensated by the equally strong y and z dipole distributions of the v y and v z field-constraint 
correlation functions. The same is true for the quadrupole distributions of the shear constraints, as 
well as for the correlation functions corresponding to the three first derivatives dfc/dxk and the 
six second derivatives d 2 fa/ 'dx^dxi. 

The predicted variance field can also be recognized when comparing the four field realiza- 
tions. They show very small differences in the neighbourhood of the peak, but further outward the 
differences become larger and ultimately are equal to the variations in any average field. 

5.3 Realizations for Gravity and Tidal Field constraints 

An important ingredient of our code is the ability to put constraints on the peculiar gravity or 
the tidal field acting on a peak. While the peaks in figure 1 do not have constraints on either the 
gravity and tidal field, we intend to give an impression of the consequences for both density and 
velocity fields of imposing such constraints by means of a sequence of four random field realiztions. 
Each of the four examples contain the same peak at the center of the box, but differ in the constraints 
on the gravity and tidal field to which the peak is subjected. The central 3uo peak is defined on 
a Gaussian scale of Ah" 1 Mpc, is spherical in shape, and has a peak curvature of V 2 /g = {x)(T2 ~ 
2.901(72. By using the same random number generator for each of the realizations we try to keep 
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the differences between the residual fields at a minimum (however, note from eq. Ill that there 
will be differences depending on which constraints are applied). 

In the first example (A, figure 5a) the central peak is not subjected to any velocity and shear 
field constraints. In the case of the second example (B, figure 5b), the same peak is constrained to 
have a peculiar velocity of 1000 km/s in the positive x-direction (note that the specified numerical 
values of these quantities are the linear extrapolations to the present epoch, a = 1, and that we 
specify the gravity and tidal field constraints in terms of peculiar velocity and shear). In the third 
realization (C, figure 5c) we constrain the shear at the peak's position, while its peculiar velocity is 
unconstrained. The off-diagonal terms of the shear tensor are zero, while a xx has a positive value 
of 50km/s/Mpc on the scale of 4.0/i~ 1 Mpc and a yy and a zz have equal and negative values. In 
the final, fourth, realization (D, figure 5d) we combine the constraints to the peculiar velocity in 
the second example and the shear at the position of the peak in the third example. 

The density and velocity field realizations for the four different constraint sets are the subject 
of figure 5. In all four cases we use a set of six panels to highlight different aspects of the fields, with 
each panel illustrating a density or velocity field in the same 5/i _1 Mpc planar section along the 
z-direction, centered halfway in the simulation box. The different contributions to the constrained 
density fields are shown in the top row panels, in combination with the corresponding velocity fields 
in the bottom row. The top left panel contains the contourmap of the mean density field, smoothed 
by a Gaussian filter with a scale of 2/i _1 Mpc (contour spacing is equal to 0.65=0.376cto(2/i -1 Mpc)). 
The corresponding mean peculiar velocity field is represented by the vector velocity field in the panel 
below. The arrows are the projections of the velocity vectors, for presentation purposes we limit 
outselves to show them at the positions of the gridpoints of a 32 3 grid. The length of each arrow 
is proportional to the magnitude of the velocity, a length of l/20th of the boxlength corresponding 
to a velocity of 1000 km/s. The corresponding full density field realization is represented by two 
panels, a density contour map of the density field (top middle panel), smoothed on the constraint 
scale of 4/i _1 Mpc, and a Zel'dovich particle distribution (top right panel). The constraints will 
heavily influence the wavevectors on a scale comparable to and larger than the scale on which they 
are imposed, while the smaller scale waves, responsible for the subclumps and other small scale 
features, are not very much affected due to their negligible correlation with the imposed constraints 
(compare eq. 38 and the listing of constraint kernels H(k) in Appendix F). The contourmap in the 
top middle panel (contour spacing 0. 275=0. 290cro(4/i~ 1 Mpc)) is therefore the best illustration of 
that part of the density field affected by the constraints. The particle distribution, on the other 
hand, shows the contribution of the small scale waves to the density field at highest possible 
resolution. The particle positions were obtained by using the Zel'dovich approximation to evolve 
an initial distribution of 64 3 particles to an expansion factor a = 0.4, approximately the time at 
which the maximum density fluctuation on the scale of 1 gridcell is equal to 10.0. An additional 
advantage of this particle distribution is that it provides a good representation of how the density 
field evolves deep into the quasi-linear regime. The velocity vector field in the bottom right panel 
is the unsmoothed full velocity field realization, and is closely related to the Zel'dovich particle 
distribution. The resolution of this velocity field representation is essentially that of one gridcell 
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Figure 5. Four different realizations of constrained random fields in the standard cold dark matter scenario (Q = i.O, h = 0.5). 
The constraints are specified on a Gaussian scale of 4h~ Mpc. In all cases there is the same Jq = 3<ro spherical peak, with 
standard curvature V 2 /g ~ 2.901o"2, at the center of the box. In (a) no further constraints are specified. In (b) the peak is 
constrained to move with a peculiar velocity of 1000 km/s towards the positive x-direction. In (c) the diagonal components of 
the traceless shear tensor are constrained to have the value o xx = 1 00 km/s/Mpc and a yy = cr zz = — 50 km/s/Mpc while the 
off-diagonal components arc all zero. In (d) the spherical peak has the combined velocity and shear constraints of (b) and (c). 
The four examples are illustrated by six panels. All show an aspect of the density or velocity field in the 5/i -1 Mpc thick central 
slice of the 100/i -1 Mpc box. Top left panel: the 2/i _1 Mpc smoothed density contourmap of the mean field /, contour spacing 
0.65. Top middle panel: the 4/i — 1 Mpc smoothed density contourmap of the constrained field realization /, contour spacing 
0.275. Top right panel: Zel'dovich particle distribution at the epoch for which the maximum density fluctuation is / = 10. on 
the scale of one gridcell. Bottom left panel: mean velocity vector map corresponding to mean density field /. The vectors are 
the projected velocity vectors in this plane. A vector with a length of l/20th of the boxsize represents a velocity of 1000 km/s. 
All velocity vector maps were determined on a 64 3 grid, but for presentation purposes only the vectors on the gridpoints of a 
32 3 subgrid are shown. Bottom middle panel: vector map of the constrained velocity field, Gaussian smoothed on a scale of 
4h.~ 1 Mpc. Bottom right panel: unsmoothed constrained velocity field vector map. 

in the 64 3 grid that was used to perform the constrained field calculations. Filtering this velocity 
field with a Gaussian function of radius AhT 1 Mpc yields the velocity field in the bottom middle 
panel, corresponding to the smoothed density field in the panel above it. 

A perfectly spherical density distribution around a maximum at the center of the box is ev- 
idently the mean density field in example A, with pure spherical infall characterizing the vector 
velocity field (left row of figure 5a). In a technical sense, recalling the discussion on figure 3, we 
can understand the spherical density field as the superposition of the spherical correlation func- 
tion (/ fa) (top middle panel) and three equally large contributions from the correlation functions 
(/ V^/g), (/ V?/g), and (/V^/g) (bottom left panel), whose main effect is to produce a slightly 
flatter peak. The spherically shaped peak can also be recognized in the center of the full field 
realization. However, the shape of the central clump becomes very irregular further outward from 
the center. A comparison with the Zel'dovich particle distribution shows that this clump consists 
of at least four separate subclumps. Note that the central peak, unlike the peak in the mean field, 
has a considerable peculiar motion in the negative y direction, and a small but nonzero shear. Both 
are introduced via the residual field. The absence of correlations between the small scale waves 
is well illustrated by the full velocity field in the bottom right panel of figure 5a, which besides 
spherical infall does not appear to display any additional features but the expected noise. 

The character of the field realization changes considerably by adding the extra constraint that 
the central spherical peak has a peculiar velocity of 1000 km/s in the x-direction (example B, 
figure 5b) . The presence of the central spherical peak can still be recognized in the mean density 
field and the full field realization. At the same time we see that the global matter distribution 
is sculpted into the dipolar pattern that induces the net gravitational acceleration corresponding 
to the required peculiar velocity. The mean velocity field in the neighbourhood of the central 
denstiy peak clearly reflects the required bulk motion. This local motion is part of a more global 
pattern in the velocity field, consisting of a convergence towards one point, 'attractor', in the right 
half and a an outflow pattern from the underdense regions in the left half. Besides this mean 
component, the full velocity field realization contains additional local features, clearly visible in 
the lower middle and right panel of figure 5b. Note that there are several local regions from which 
matter is streaming away, some of these local density depressions are not even underdense (note 
e.g. the saddlepoint around [x,y] = [70.0, 50.0]h~ l Mpc). Also remark the fact that the central 
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peak is more compact than in example A, mainly due to the very steep density falloff of the peak on 
the side where it is lying on the boundary of the underdense region. This pattern finds its origin in 
the extra superposition of the dipolar pattern characteristic of the correlation function (/ vq,x) (see 
figure 3). Equally striking are the consequences of imposing extra constraints, in example C, on the 
tidal field and/or corresponding shear at the peak position (figure 5c). The constraints induce the 
expected global quadrupolar mass distribution in the mean density field, superimposed on the local 
spherical peak density distribution. The band of matter parallel to the x-axis, visible in both the 
mean and final density field, induces the dilational shearing motion along the x-direction and the 
compressional shear along the other two directions, in collaboration with the underdense regions 
below and on top of it. The presence at the peak position of the positive a xx component, along with 
the negative a yy component of half its magnitude, is most strikingly visible in the mean velocity 
field. In the full field realization we can also recognize the presence of other components than the 
quadrupolar one. The central high-density ridge is littered with numerous small scale peaks of 
different sizes (see e.g the Zel'dovich particle distribution) while a clear dipolar component can also 
be discerned in the density distribution. High-density regions are concentrated in the lower half 
of the box, inducing the sizable peculiar motion of the peak towards the negative y-direction that 
can be seen in the velocity field realizations in the lower middle and right panels. Finally, figure 5d 
shows how the combination of the constraints on the peculiar velocity and the shear in example D 
work out. The corresponding mean density field clearly contains both a dipolar and a quadrupolar 
component, both of which are also conspicuously present in the full density field realization (also 
compare with the Zel'dovich particle distribution). In both the mean velocity field and the full 
velocity field realization we can recognize the specified peculiar velocity and shear at the position 
of the central peak. The particle distribution shows that the clumps on the right hand side of the 
center are more massive than the ones in figure 5c. The agglomerate of these clumps conspires to 
form a big attractor, easily recognizable, that induces the large peculiar motion of the peak. 

6. Summary and Discussion 

In this paper we have developed a formalism to set up cosmological initial Gaussian random density 
and velocity fields that can contain one or more peaks or dips, with the intention to generate 
appropriate initial conditions for cosmological iV-body simulations that focus on the evolution of 
the progenitors of the present-day galaxies and clusters and their environment. The method is suited 
for fields with any arbitrary power spectrum P(k). Central objective of our algorithm is the ability 
to sculpt the local and global matter distribution in a sufficiently large volume such that certain 
physical characteristics of the density and velocity field in the immediate neighbourhood of the 
primordial peaks have a priori specified values. The generation of these constrained density fields is 
an application and elaboration of the the Hoffman & Ribak (1991) prescription. They showed that 
there is a simple and elegant solution to achieve this if the constraints are linear functionals of the 
field. We have presented the implementation of our method following a comprehensive discussion 
of the fundamentals underlying their method. 

A maximum of 21 characteristics is used to specify the density and velocity at and around the 
position of the peak. They can be divided into three groups: 
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[1] The scale and position of the peak. We identify a peak as a local maximum in the density 
field that has been smoothed by a Gaussian filter function with a characteristic scale Rg, although 
the formalism can be very easily extended to other filter functions. The peak may be positioned at 
an arbitrary position within the simulation box. 

[2] The local density field. In total 10 constraints are needed to fully specify the density field 
in the immediate vicinity of the peak. The first one concerns the height of the peak. In addition, 
three constraints are needed to assure that the three first derivatives of the smooth density field 
vanish at is summit. Finally, the six second order derivatives of the smooth density field are set by 
specifying the compactness V 2 /g> the axis ratios and the orientation of the peak. 

[3] The local gravitational field. The specification of the gravitational field around the peak 
introduces 8 additional constraints: the three components of the smoothed peculiar gravity at the 
location of the peak and the five independent components of the traceless tidal tensor. The resulting 
density field is sculpted in such a way that it induces the desired amount of net gravitational and 
tidal forces. We usually restrict ourselves to the growing mode component of the density field. In 
the linear clustering regime the peculiar gravity and tidal field are therefore directly proportional to 
the peculiar velocity and the shear, so that we commonly use the latter to specify the gravitational 
field constraints. 

It may be worthwhile to point out that in a linear density fluctuation field several of the above 
quantities are correlated. For example, we find that there is a strong correlation between the tidal 
field tensor and the mass tensor, expressing itself in the tendency of the tidal field to align itself 
along the principal axes of the mass tensor. 

The constraints that we consider here are linear functionals of the density fluctuation field /, 
and therefore can be written as convolutions of the field with a specific function. Consequently, 
it is most convenient to perform the relevant calculations in Fourier space. The generation of 
a constrained field realization basically consists of the sum of an arbitrary field realization with 
the convolution of the power spectrum with a function that is the weighted sum of the different 
constraint kernels, the weights depending on the specified values of the constraints and the values 
of the constraints for the unconstrained field (see eq. 26). The expressions for these constraint 
kernels are derived from the particular constraints to which they are related. In Appendix F we 
list the kernels used in our code. 

The Hoffman-Ribak algorithm that we have described here is considerably faster and more 
generally applicable than the original Bertschinger (1987) algorithm. Its superior speed is due 
to the direct and simple way of sampling the residual field, rendering an iterative "simulated 
annealing" technique superfluous. Moreover, because it is a direct method it has the additional 
advantage of superior accuracy. Extensive testing of constrained field realizations showed that 
the implementation is very precise, leading to accuracies in the order of 0.01% for the imposed 
quantities. In the computer implementation of our code the constrained field is evaluated on a 
periodic three-dimensional lattice. This has the advantage of being able to perform the Fourier 
transforms by means of a Fast Fourier Transform, with the advantage of being considerably faster 
than methods based on a direct Fourier transform. A disadvantage of the FFT is that they have 
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a rather weak sampling at low k, while direct Fourier transforms enable a far better sampling in 
that range. In their multipole constrained field method Bond and Meyers (1993) therefore resort 
to direct Fourier transforms, resulting in an excellent sampling at low and intermediate k. 

In addition to the fact that the Hoffman-Ribak method provides us with a fast, efficient and 
accurate method to generate constrained random fields it has two other important advantages. The 
first one is that the implementation of a large variety of constraints is relatively straightforward 
through the convolution integrals in Fourier space. Secondly, unlike most other efficient algorithms 
it is equally suitable and efficient for local and non-local constraints. Although the illustrations of 
the peak constraints in section 5 were mainly local in character, centered on one peak, the developed 
formalism allows the generation of numerous peaks and dips at different positions (see figure 1). 

In our application to peaks we followed the philosophy that each of the constraints corresponds 
to a different physical quantity. Another class of possible applications of the Hoffman-Ribak pro- 
cedure is the reconstruction of (linear) density fields from the measurement of the same physical 
quantity at several different positions inside a certain volume. A nice illustration of this is the work 
by Ganon & Hoffman (1993), who reconstructed the the density field in the "local" universe from 
the observed velocity field sampled at 181 different positions within a sphere of 40/i~ 1 Mpc around 
us, assuming that it is a realization of a standard cold dark matter field. They showed that the 
method recovers the main features of POTENT's density field (Dekel, Bertschinger & Faber 1990), 
in particular the Great Attractor region. The interesting feature of this reconstruction application 
is that it creates high-resolution fields subject to the low-resolution data, for the given underlying 
model. It therefore offers the charming and interesting opportunity to set up initial conditions for 
A-body simulations from observations of the local Universe, so that the nonlinear evolution of our 
"local" Universe in a particular cosmological scenario can be studied. A related and promising 
application would be the construction of high-resolution microwave background maps from the the 
large-scale anisotropies measured by COBE (Bunn et al. 1994). 

This class of constraint problems, where the constraints consist of the value of the same phys- 
ical quantity Y>(r) a ^ many different positions, offers the advantage that for every constraint the 
constraint-field correlation function £j(r) = (ip(ri)f(r)) = T(r — r,) (see eq. 24) can be evaluated 
from the same general correlation function T(x). The same is true for the constraint-contraint 
correlation function In particular, this will be a great advantage if the constraint values are 
imposed at equally spaced points on a grid. This is the approach followed by Ganon & Hoffman 
(1993), who determined the constraint values for the velocity potential on a grid by spatial in- 
terpolation from observed values of the peculiar velocity. The computation of the required values 
of £j(x) and the inverse constraint-constraint correlation matrix ^ 1 can then be simply accom- 
plished by two FFTs. This can be easily seen from the following. Because the quantity V( x ) is a 
linear functional of the density field /(x), its Fourier transform ip(\t) is a product of the Fourier 
transform /(k) of the field /(x) with a kernel function h(k), V>(k) = /i(k)/(k). Examples of such 
fields V( x ) are the gravitational potential, the peculiar velocity in the linear regime, or the tem- 
perature variations in the cosmic background radiation field. After evaluating the corresponding 
expressions for h(k) at wavenumbers k p (compare the kernel functions listed in Appendix F), the 
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values of £i( x j) = (^(xj)/(xj)) and £y = (V'(xj)V ; (x : , )) can be found from £i(xj) = T(xj — Xj) and 
&j = *(xj - Xj), where 



T(x) = l^(k p )%)e^ 

P= ° • (H3) 

*( X ) = ^E Mk P ) 2 P(Me-^' x 

p=0 

In fact, the inverse matrix of £jj can be found directly and very simply from = 0(x/j — xj), 
where 9(x) is the inverse of *(x), i.e. ^(x^ — Xj)0(xj — x^) = 5^, and therefore given by the 
Fourier sum 

1 JV " 1 1 

9(x) = —Y = e~^' x . (114) 

The computation of the discrete Fourier sums T(x) and B(x) is accomplished by a FFT, so that 
the computational cost is only O(NlogN). Note that because of the periodic boundary conditions 
intrinsic to the FFT each coordinate of Xj can only attain half of the values along each axis, so 
that in total only | of the computational box is used for the field reconstruction. Finally, the 
independent Fourier components of the unconstrained field /(x) are generated. The subsequent 
computation of /(x) itself demands one FFT, and the computation of the corresponding constraint 
values Cj involves another FFT (compare eq. 38) . Combining all these results in the final evaluation 
of the constrained field according to equation (35) consists of the computation of the double product 
^(x)^ 1 (cj — Cj) for every point Xj, making it an 0(N 3 ) formalism. However, unlike the formalism 
developed in section 3, this procedure does not involve a very costly matrix inversion of £y , implying 
it to be far more efficient and the method of choice for this particular class of applications. On the 
other hand, when each of the M constraint quantities have a different character, concern different 
scales, arbitrary non-grid positions, or different filters, this procedure cannot be straightforwardly 
applied. In those cases a formalism similar to the one presented in this paper is automatically 
implied. 

As a final note we should issue a cautionary remark on the practical implementation of our 
constrained random field code. The initial density fields are set up in a box with periodic boundary 
conditions. This means that the mean density of the box is exactly equal to the mean density of the 
Universe. The structure generated within the box is therefore not entirely typical, since overdense 
regions must necessarily be surrounded by low-density regions. This need not be true in general, 
from the theory of Gaussian random field we know that peaks tend to cluster. The simulation box 
should therefore not be taken too small, the resulting structure might be very atypical. Evidently, 
this conflicts with the demand to make the box as small as possible to achieve the highest possible 
resolution. The chosen box size should therefore be a compromise between these two. 

In summary, we can conclude that the Hoffman-Ribak method provides a powerful and elegant 
tool to study the formation and evolution of specific cosmological objects in great detail under 
ideal conditions. The tools developed in this paper should essentially be regarded to constitute a 
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laboratory equipment set*. They allow us to set up very specific conditions for the objects under 
study. A sequence of experiments based on a range of different circumstances will subsequently yield 
a maximum of insight into the systematic dependence of structure formation on specific physical 
quantities. By concentrating on one specific application, peaks in the density field, we hope to 
have provided a recipy for constructing similar applications and extensions for different quantities 
in fields of a possibly different character. A straightforward extension of our formalism will for 
example be to consider peaks in the gravitational potential field instead of in the density field. 
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Appendix A. The intersection of a sphere and a polygon 

In section 2.2, equation 19, we saw that imposing the set of constraints T = |Cj[/;Xj] = cf, i = 
1, . . . , M} is equivalent to a change of the action S[f] into 

2S[f] = j J /(xi)tf(xi - x 2 )/(x 2 ) dx 1( ix 2 - C^C, (Al) 

with £ij the (ij) th element of the matrix Q = (C*C). 

The constraint Cj(xj) (in this appendix we will use the simplifying notation Cj(x) for 
Cj[/;x] ) can be written as a convolution with a Dirac delta function 5d(x), 

Cj(xj) = J dx 2 d D (yi2)Ci(xi -x 2 ) = / ^x 2 £( x i)^( x i - x 2 )Q(xj - x 2 ) , (A2) 

where we have used the fact that K(x) is the functional inverse of the correlation function £(x) 
(eq. 7, section 2.1). By using the convolution theorem we can express this double convolution 
integral in Fourier space as 



Cite) = / 7^3 Q(k)P(k)^(k) e-* k - , (A3) 



where C{ (k) is the Fourier transform of Cj (x) , 

r rl\ 

Q( x ) = / ^a(k)e- h , (A4) 

and -P(k) = P(&), the spectral density, and K(k) the Fourier transforms of £(x) and K(x) respec- 
tively (eq. B5). The formal definition for the spectral density P(k) is (Bertschinger 1992) 

(27r) 3 P(fci) <5 D (k! - k 2 ) = (/(ki)/*(k 2 )) , (A5) 

with <5o(ki — k 2 ) the Dirac delta function. In an analogous fashion a function Pj(k) can be 
introduced, 

(2^) 3 p(kx) ^(ki - k 2 ) = (Q(k0r (k 2 )) , (A6) 
from which we obtain, in combination with equation (A5), a relation between /(k) and Cj(k), 

(/(kQ/^k,)) _ (cUkQnkg]) rr^-^ffki (An 

By subsequently inserting this relation in the Fourier integral of equation (A3), and using defini- 
tion (A6), we get 
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Ci(*i)= J -|-3/(k)P,(k)^(k)e^ 



(27T) 

= // (^ja ( ^)3 /(ki)^(k 1 )(a(k 1 )r(k 2 ))e- k1 ^ , (A8) 
= y J /(xi)K(xi - x 2 )^(x 2 ) dxi dx 2 

where the function £j(x) is the field-constraint correlation function and the Fourier transform of 
£(k), 

fc(x) = (/(x)Q(x,)) = / ^3 ^(k)e- k -(—) . (A9) 

Since the field /(x) also obeys the constraints Cj = Cj the expression Cj(xj) £^ 1 Cj(x ? -) in equa- 
tion (Al) can be replaced by 

C< f^Ci = / / /(xi)^(xi - x 2 )/(x 2 ) d Xl dx 2 , (A10) 
where we have defined the field /(x) by 

7(x) = £(x)^V (All) 
The constrained action in equation (Al) can therefore be written as 

2S[f}= J J /(x 1 )K(x 1 -x 2 ){/(x 2 )-/(x 2 )}dx 1 dx 2 

= J J {/(xi) - 7( Xl )}K( Xl - x 2 ){/(x 2 ) - 7(x 2 )} dx lC /x 2 + (A12) 
1 1 7(x 1 )K(x 1 -x 2 ){/(x 2 )-7(x 2 )}dx 1 dx 2 



4-2 • 



It can be easily shown that the second term on the right hand side of equation (A12) is equal to 
zero because 



J J /(xi)K(xi - x 2 )/(x 2 ) dxidx 2 = j J /(xi)i^(xi - x 2 )/(x 2 ) dxidx 2 = (H^Cj . (A13) 

This relation follows directly from equation (A10) for the second integral, while for the first integral 
it follows from the fact that 



J J /(xi)-ff(xi - x 2 )/(x 2 )cfocidx 2 = Zi/tuCjCi J J ^(xi)if(xi - x 2 )£ fc (x 2 ) dxidx 2 



= Zij^utkiCjCt = ci i l3 x Cj 



(A14) 



where we have used the substitution of the Fourier expression of &(x) (eq. A9) to evaluate the 
integral, 
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J J Ci( x i)^( x i - x 2 )Cfc(x 2 ) dxidx 2 



^3 P,(k 1 )K(k 2 )P fe *(k 2 ) (2vr) 3 fe(k 2 - kx) e-^e^ (A15) 

= //^^(C fe *(k 2 )a(k 1 ))e- ik ^e^ = (C fe (x fe )a(x,)) = . 

The transition from the 2 nd to 3 rd line in equation (A15) has been made by combining equation (A6) 
and (A7), 

(2vr) 3 A(k 1 )fe(k 1 - k 2 ) = |M (C^C^)) , (A16) 

and the fact that P{k) = l/K(k) (see eq. B6, app. B). 

By defining the "residual field" -F(x) = /(x) — /(x) we can therefore conclude from equa- 
tion (A12) that the constrained action S[f] = S[F] can be written as 

2S[F] = J Cbd J dx 2 F(x 1 )if(x 1 - x 2 )F(x 2 ) , (A17) 
which is the expression needed in Section 2.2. 

Appendix B: Diagonalisation of the action S[F] 

In this appendix we will rewrite the action S[F] (eq. 20), 

S[F] = lfdxif dx 2 F*( Xl )K( Xl - x 2 )F(x 2 ). (Bl) 
in terms of the Fourier transform F(k) of the fluctuation field F(x), 

r d\i 

F(X) = J (27f F(k)e " lk ' X - (B2) 

The kernel K(x) in equation (Bl) is the functional inverse of the correlation function (eq. 7), 

J dx^( Xl - x)£(x - x 2 ) = 5 D (xi - x 2 ) . (B3) 
By virtue of the convolution theorem this equation is equivalent to 

J ^ K(k)P(t)e fc ^) = fo( Xl - x 2 ), (B4) 

where K(k) and P(k) are the Fourier transform of K(x) and £(x) , 

r dk r r/k 

The identification of the left part of equation (B4) with the Fourier integral expression of the Dirac 
delta function implies that K(k) = 1/P(k) . Consequently, 
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f dk 1 



(B6) 



Likewise, the insertion of K(k) = 1/P(k) into the double convolution of equation (Bl), 




(B7) 



yields the Fourier expression for S[F], 



S[f] 



I 



dk |F(k)| 2 



(B8) 



(2vr) 3 2P(k) 



which is equation 29 in section 3. 

Appendix C: A heuristic proof that V [F\T] is independent of q. 

A field /(x) can be viewed as an iV-dimensional vector (/i, . . . , /jv) in iV-dimensional "field" space, 
with N — > oo. The fields /(x) that obey the set of M constraints T = {Ci[f; Xj] = q; i = 1, . . . , M} 
define an (iV— M)-dimensional hypersurface in this iV-dimensional space. For reasons of convenience 
this hypersurface will also be denoted as T. The only restriction that we impose on the constraints 
Ci is that they are linear, 



Each of the hypersurfaces T contain a special point /(x), the mean of the fields satisfying the 
constraints T, 



where £i(x) is the cross-correlation between the field and the i th constraint Cj[/;x], and £jj the 
correlation between the i th and j th constraints, Ci[f] and Cj[f] (notice that in this notation we 
stress the functional character of the constraints). Both £jj and £«(x) are defined in equation (24) 
(Section 2.2). Each of the fields /(x) in T have a corresponding residual field -F(x), defined as the 
difference between the field /(x) and the mean field /(x) of T, 



Imagine two arbitrarily chosen constraint hypersurfaces, the first one corresponding to the con- 
straint set Ti = {Cj[/;xj] = Cj ; i; i = 1,...,M} and the other one to the set T2 = {Cj[/;xj] = 
Ci,2'i i = 1, • • • , M}. The mean fields of the sets F\ and T2 are fi and fo. Consider the translation 
of an arbitrary field /i(x) G T± by a field T2 j i(x) into a field /r(x), 



C i [/i + / 2 ;x]=C i [/i;x] + C i [/ 2 ;x]. 



(CI) 



/(x) = (/(x)|r)=^(x)^ 1 Cj 



(C2) 



F(x) = /(x)-/(x). 



(C3) 



where the translation T2 5 i(x) is defined by 
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(C4) 



T 2)1 (x) = / 2 (x) - /!(x) = &(x) &( Cj , 2 - c hl ) . (C5) 

This definition of T 2t \ immediately implies that the mean field /i(x) of T± is transformed into the 
mean field / 2 (x) of r 2 . From equations (C4) and (C5) and the linearity of the constraints Q we 
can infer that 

C i [f T ]=C i [f 1 ] + C i [T 2tl ] 

= C i [h\ + C i [f 2 \-C i [h\ (C6) 

= Ci,l + Ci,2 - Cj,l = Cj,2 , 

The field /t(x) therefore obeys the constraint set T 2 . This is true regardless of the field /i(x) G Ti. 
Moreover, the inverse translation — T 2j i transforms the resulting field /2(x) back into /i(x). The 
two hypersurfaces T\ and T 2 are therefore linked by a one-to-one mapping, so that 

nfi|ri] = n/2|r 2 ], (C7) 

where "P[/i |Ti] is the probability of having a specific field /i(x) under the condition that they 
satisfy the constraints Ti, and / 2 (x) is the field in the hypersurface r 2 that is linked to /i(x) by 
the translation T 2i i(x) (eq. C4). The conditional probabilities for the corresponding residual fields 
i*i(x) = (/i(x) — /i(x)) and F 2 (x) = (/ 2 (x) — / 2 (x)) can be inferred from equation (C7), 

P[Fi|ri] = rihlFi] = r[f 2 \r 2 ] = v[F 2 \r 2 ] . (cs) 

Finally, consider the transformation of the residual field -Fi(x) under the translation T 2i i, 

F 1 (yi) = f 1 ( X )-f 1 (yi) 

= (A(x) + T 2il (x)) - (7x(x) +T 2il (x)) (C9) 
= / 2 (x)-/ 2 (x) = F 2 (x). 

In other words, the residual field -F(x) is invariant under the translation T 2j i, i.e. F± = F = F 2 , 
which in combination with equation (C8) implies that 

p[F|ri] = pfFxirx] = v[F 2 \r 2 ] = p[f\t 2 ] . (cio) 

This is the result that we intended to prove. 

Appendix D: The variance (F 2 (x)|T) of the residual field F(x). 

A derivation will be given for the expression for the variance (i ?2 (x)|r) of the residual field belonging 
to the constraint set T. The residual field F(x) is the difference between a field /(x) obeying the 
constraint set T and the mean /(x) = (F(x)|T) of all these fields, 



7(x) = ^(x)^ 1 c J - 
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(Dl) 



The crucial observation that V [F\F], the probability of having a residual field -F(x) satisfying a 
particular set of constraints T, is independent of the numerical value Cj of the constraints T, 



V \F\Y{\ = V [F|r 2 ] , (D2) 

implies that 

<F 2 (x)|r> = <F 2 (x)), (D3) 

where (F 2 ) is the variance in all possible realizations of the field, and (F 2 \F) the variance for the 
ones that obey the constraint set T. From equation (D3) we find 

<F 2 (x)> = (V[T] <F 2 (x)|r> = [V[T] <(/(x) -7(x)) 2 |r> 

(D4) 

= /p[r] |(/ 2 (x)|r)-(/(x)|r) 2 }, 

where V [r] is the integrated probability of all realizations that obey the constraint set V. Evaluation 
of the first part of the integral in (D4) yields 

/ v [r] </ 2 (x)|r> =\v\v\\v [/(x)|r] / 2 (x) = jv [r] v [/(x)|r] / 2 (x) 

= jr[f(x)] / 2 (x) = (/ 2 (x)) = a 2 , 

where a 2 , is the general variance of the density field fluctuations. In the derivation of (D5) we have 
used the fact that P[/|r] is the product of the probability V [r] with the conditional probability of 
having the field /(x) under the condition that it obeys V, V [f\T] (equation 13, section 2.2). 

To evaluate the second part of the integral we use the expression for the mean field in 
equation (Dl), 

JV[T] (/(x)|r) 2 = JV[T] ti(x)&c jCl Z-%(x) 

= fcCx)^ 1 {/ V (T) Cj^Qfr)} Cr/6(x) ( D6 ) 

= ^(x)^ 1 (c j Q)er fe 1 6(x) = ZiW&ZjitikZkW = tiW&tkW. 

By inserting equations (D5) and (D6) into equation (D4) and using equation (D3) we find 

(F 2 (x)|r)=a 2 -e i (x)^(x), (D7) 

which is the intended expression. 

Appendix E: Shape and orientation of a peak in a random field. 

The second order Taylor expansion of a density field around a peak or dip at position x^ in a 
density field /(x) is given by equation (51), which we repeat here for convenience, 
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1 3 d 2 fa 

/g(x) = /c(x d ) + ^ S a a . ( x rf) ( x * - ~ x d,j) • (El) 

This quadratic equation can be written in its canonical form by transforming to the coordinate 
system x' = {x^, x' 2 , whose axes are aligned along the eigenvectors of the matrix VjVj/g- If 
the eigenvalues of — VjVj/ 9 are Ai, A2 and A3, equation (El) becomes 

/g(x') = / o (0) - i £ A*x? , (E2) 
z i=i 

where we have chosen the origin of x' to coincide with the position of the peak or dip. In the case 
of a peak the Aj have a negative value, for a dip they have a positive value. From equation (E2) 
we see that the isodensity surface fa = F is a triaxial ellipsoid whose principal axes are oriented 
along the coordinate axes, with semiaxes given by 



In equation (E3) the central height /g( x (z) of the overdensity is expressed in units of ao(Ra), i.e. 
fc{*d) = v c a (R G ). 

From equation (E3) and the fact that the shape of a triaxial ellipsoid is fully specified by its 
two axis ratios a\i = (01/02) and 013 = (01/03) we can infer that constraints on the shape of the 
overdensity result in constraints on the ratio of Aj's, 

_ 2 f ^3 



xj IatJ=°- < E4) 

The actual magnitude of the Aj's depends on the steepness of the density profile around the peak. 
This steepness is specified by the Laplacian V 2 /g, as can be observed from the expansion of the 
density profile equation (E2) in spherical coordinates (x,6,ip), 

2 

/ G (x) = / G (xd) + V 2 /G(x rf ) y {1 + A(9, <p)} . (E5) 

A(6, ip) is a function of the direction (9, (p) that describes the asphericity of the peak via its depen- 
dence on the parameters Ai, A2 and A3. In deriving equation (E5) we used the relation between 
the A, and V 2 fc 

3 

E^ = - V2 /G(x rf ), (E6) 
i=i 

which can be obtained by double differentiation of equation (E2). Usually V 2 /g is expressed in 
units of U2{Rg) = (V 2 /g V 2 /g) 1//2 , i.e. V 2 /g(-Rg) = — Xd^iRc)- The expression for A 1 is obtained 
by combination of equations (E4) and (E6), 

Ai = - X J^M_ . (E7 ) 
(1 + o 2 2 + o 2 3 ) 
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Once the value of Ai has been determined, the values of A2 and A3 are obtained by multiplication 
of Ai by a\ 2 and af 3 respectively (equation E4). 



The orientation of the peak with respect to the general coordinate system is described by the three 
Euler angles a, (3 and tp (see Goldstein 1980). Here (3 is the angle between the smallest axis of 
the ellipsoid and the z-coordinate axis, a the angle between the line of nodes and the x-coordinate 
axis, and tp the angle between the largest axis of the ellipsoid and the line of nodes. The line of 
nodes is the intersection of the xy-plane and the plane defined by the largest and second largest 
axis of the ellipsoid. The transformation matrix Aij (equation 60, Sect. 4.2) is obtained from this 
definition of the Euler angles, 



cos a cos tp — cos (3 sin a sin tp sin a cos tp + cos (3 cos a sin tp sin (3 sin tp 
A = I — cos a sin tp — cos (3 sin a cos ip — sin a sin tp + cos (3 cos a cos ip — sin (3 cos tp ] . (E8) 
sin (3 sin a — sin (3 cos a cos (3 

This matrix describes the transformation from the coordinate system x' defined by the principal 
axes of the ellipsoid to the general coordinate system x, 

3 

x'i = ^2 Aijixj - x d:j ), i = l,..., 3, (E9) 
i=i 

with Xd the position of the centre of the peak. Thus, xf transforms as 

3 3 

X 'i = Y Y A ij A 'ik{Xj ~ X d j)(x k ~ X d:k ) . (E10) 
j=l k=l 

By inserting this transformation into the expression for the density profile (eq. E2) we obtain the 
following quadratic equation for the density profile in the general coordinate system x, 

1 3 f 3 1 

/g(x) = fa(xd) - j E E X i A ij A ik \ (xj - x dJ )(x k - x d>k ) . (Ell) 

j,k=i U=i J 

Because equation (Ell) is equivalent to equation (El) we obtain the following relationship between 
the second derivatives of fa and the orientation, shape and steepness of the dip or peak in the field 

(> ' J(; -J2Xk A ki A kj, i, j = 1,2,3. (E12) 



This is the expression that we use in section 4.2. 
Appendix F: Peak constraint kernels and values. 

Here we present the explicit expressions for the 18 peak constraints and the corresponding kernels 
Hl(k), defined in equation (37), to give an overview and summary of the results in this paper. 
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The filter function W(k) is taken to be the one corresponding to a Gaussian filter function with 
smoothing length Rq, 



W(k) 



-k 2 R 2 j2 



(Fl) 



The peak constraints are presented in 5 groups. The first group consists of the peak height constraint 
/G*( x d)) the second one of the three constraints on the first derivative of the field, Vf G (xd), and 
the third one of the second derivatives VjVj/G(xd)- In addition, the fourth group contains the 
constraints on the peculiar velocity of the peak, vg(x^), while the fifth group corresponds to the 
constraints on the five components of the shear, o"G,ij(xrf), 



/c(xd) = v(Tq{Rg) 



d 2 f G 

d 2 j\ 



¥ (xd) = -ELi ^A kl A kl 



9 2 .fG 
dx\dx2 

d 2 f G 



^T'(xd) = - Efc=l AfcA fc3 A fc3 

(Xd) = - Efc=l ^kAklAk2 



^(k) = e -fc 2 «y2 e i k .x d 



# 3 (k) = ifc 2 e - fe2R G/ 2 e ik - x " 
fl- 4 (k) = ifc 3 e^ R H 2 e ik ' Xd 



ff 5 (k) = _fe 2 e - fc2 R G / 2 e ^ 

H 6 (k) = -kle-^^Gl 2 ^^ 

H 7 (k) = -k 2 e- k2R y 2 e ik ^ 

H 8 (k) = -k 1 k 2 e- k2R y 2 e i ^ 
H 9 {k) = -jfe 1 fc 3 e- fc2j %/ 2 e ik - x <' 
fl"io(k) = -fe 2 fe 3 e- fc2fl G/ 2 e ik - x d 



9G,i(x<i) = 


: 9l&g,pk(R G ) 






= §Off 2 


ffG,2(xrf) = 


'- g2 a g,pk{R G ) 






= §Off 2 


ffG,3(x rf ) = 


'- 9^g,pk{R G ) 




Hi 3 (k) 


= §Off 2 


E G ,n(*d) 


= ~ea E (R G ) EL 


=1 Q(^)TkiT kl 


^14 (k) 


= I^IH 2 


^G,22(x d ) 


= ~ea E {R G ) EL 


-_i Q(ro)7fc 2 Tfc2 


#i 5 (k) 


= I^IH 2 


^G,12( x d) 


= ~ea E (R G ) EL 


= i Q(ro)r fel r fe2 


&ie(k) 


= 3qh 2 


^G,13(x d ) 


= €(Te{Rg) EL 


= i Q(ro)T fc iT fc3 


^ir(k) 


= I^IH 2 


#G,23(Xd) 


= €<J E {Rg) Efe= 


=1 Q(ro)7fc2Tfe 3 


ffis(k) 


= %nH 2 



2 ikl e -k 2 R 2 Q /2 e ik-x d 

ifo e -k 2 R 2 G /2 e ik-x d 

e -k 2 R 2 G /2 e ik-x d 



F 
F 



-k 2 R 2 G /2 ik-x d 



1) e -fc 2 R G /2 e ik-x d 



3^2 (fclM e -fc^y 2 glk . Xd 



fclfc 3 

~F^ 



-k 2 R? G /2 e ik-x d 
-k 2 R 2 Q /2 e ik-x d 
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